rnaseq/step1_fastqc/snippet.sh

@@ -1,12 +1,5 @@
-#!/bin/bash
+# The command to do a FastQC on a fastq file is
-R1="sample1_R1.fastq.gz"
+file="file_to_analyse.fq.gz"
-R2="sample1_R2.fastq.gz"
+fastqc_out="./path/to/fastqc/output/dir"
 
-PROJECT_DIRECTORY="/groups/umcg-griac/tmp01/rawdata/$(whoami)/rnaseq"
+fastqc -o "$fastqc_out" "$file"
-FASTQC_OUT="${PROJECT_DIRECTORY}/step1/"
-mkdir -p "${FASTQC_OUT}"
-
-# Run FastQC on paired-end data.
-fastqc \
-    -o "${FASTQC_OUT}" \
-    "${R1}" "${R2}"

rnaseq/step2_trim/snippet.sh

@@ -1,12 +1,9 @@
 #!/bin/bash
-#
+
-# Reference: http://www.usadellab.org/cms/?page=trimmomatic
+# reference: http://www.usadellab.org/cms/?page=trimmomatic
 
 
 module load Trimmomatic
-PROJECT_DIRECTORY="/groups/umcg-griac/tmp01/rawdata/$(whoami)/rnaseq"
-FASTQ_OUT="${PROJECT_DIRECTORY}/step2/"
-mkdir -p "${FASTQ_OUT}"
 
 
 # Adapters can be found at
@@ -14,26 +11,16 @@ mkdir -p "${FASTQ_OUT}"
 # But should be verified with FastQC, or in another way.
 
 
-# Trimmomatic example Paired end data.
+# (Example) Paired end
-#
-# Flags:
-# - ILLUMINACLIP: Cut adapter and other illumina-specific sequences from the
-#                 read.
-# - SLIDINGWINDOW: Perform a sliding window trimming, cutting once the average
-#                  quality within the window falls below a threshold.
-# - LEADING: Cut bases off the start of a read, if below a threshold quality.
-# - TRAILING: Cut bases off the end of a read, if below a threshold quality.
-# - HEADCROP: Cut the specified number of bases from the start of the read.
-# - MINLEN: Drop the read if it is below a specified length.
 java -jar $EBROOTTRIMMOMATIC/trimmomatic.jar PE \
   -phred33 \
-  sample1_R1.fastq.gz \
+  input_forward.fq.gz \
-  sample1_R2.fastq.gz \
+  input_reverse.fq.gz \
-  "${FASTQ_OUT}/sample1_R1_paired.fastq.gz" \
+  output_forward_paired.fq.gz \
-  "${FASTQ_OUT}/sample1_R1_unpaired.fastq.gz" \
+  output_forward_unpaired.fq.gz \
-  "${FASTQ_OUT}/sample1_R2_paired.fastq.gz" \
+  output_reverse_paired.fq.gz \
-  "${FASTQ_OUT}/sample1_R2_unpaired.fastq.gz" \
+  output_reverse_unpaired.fq.gz \
-  ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
+  ILLUMINACLIP: TruSeq3-PE.fa:2:30:10 \
   LEADING:3 \
   TRAILING:3 \
   SLIDINGWINDOW:4:25 \
@@ -41,11 +28,11 @@ java -jar $EBROOTTRIMMOMATIC/trimmomatic.jar PE \
   MINLEN:50
 
 
-# Example single end data.
+# (Example) Single end
 java -jar $EBROOTTRIMMOMATIC/trimmomatic.jar SE \
   -phred33 \
-  sample1.fastq.gz \
+  input.fq.gz \
-  output.fastq.gz \
+  output.fq.gz \
   ILLUMINACLIP:TruSeq3-SE:2:30:10 \
   LEADING:3 \
   TRAILING:3 \

rnaseq/step3_align/snippet.sh

@@ -2,42 +2,32 @@
 #
 # Align reads against reference genome.
 
-REFERENCE_DATA="/groups/umcg-griac/prm03/rawdata/reference/genome"
+STORAGE="/groups/umcg-griac/tmp04/rawdata/$(whoami)/step3"
-
+# Store genome index in this location:.
-PROJECT_DIRECTORY="/groups/umcg-griac/tmp01/rawdata/$(whoami)/rnaseq"
+GENOME_INDEX="${STORAGE}/genome_index"
-# Store the generated `sample1_Aligned.sortedByCoord.out.bam` in this dir.
+mkdir -p "${GENOME_INDEX}"
-ALIGNMENT_OUTPUT="${PROJECT_DIRECTORY}/step3/alignment"
+# Store the generated `Aligned.sortedByCoord.out.bam` in this dir.
+ALIGNMENT_OUTPUT="${STORAGE}/alignment"
 mkdir -p "${ALIGNMENT_OUTPUT}"
 
-# 1) Generate genome index (optional).
+# 1) Generate genome index.
-#
-# Depending on your read size, reference genome and annotation, you may need to
-# generate a new genome index. In most cases, this is not necessary and you can
-# directly use the pre-build genome index from the cluster:
-#
-# GENOME_INDEX="${REFERENCE_DATA}/index_GRCh38_gtf100_overhang100"
-#
-# and ignore the first STAR command below.
 #
 # N.B.:
-# - We're assuming a read size of 100 bp (--sjdbOverhang 100). An alternative
+# - We're assuming a read size of 100 bp (--sjdbOverhang 100). Refer back to the
-#   cut-off is 150, for low-input methods. In general, refer back to the
 #   previous quality control steps if you are unsure about the size. In case of
 #   reads of varying length, the ideal value is max(ReadLength)-1.
-# - If you're using gzip compressed reference data, i.e., .gtf.gz and fa.gz,
+# - We're using gzip compressed reference data (--readFilesCommand zcat), i.e.,
-#   pass the `--readFilesCommand zcat` flag.
+#   .gtf.gz and fa.gz. If not, you can remove the `zcat` flag.
 
-# Store created genome index in this location:.
+# Storage location reference data (in this case on calculon).
-GENOME_INDEX="${PROJECT_DIRECTORY}/step3/genome_index"
+REFERENCE_DATA="/groups/umcg-griac/prm02/rawdata/reference/genome"
-mkdir -p "${GENOME_INDEX}"
+GTF_FILE="${REFERENCE_DATA}/Homo_sapiens.GRCh38.100.gtf.gz"
-
+FASTA_FILE="${REFERENCE_DATA}/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz"
-# Storage location reference data on Gearshift.
-GTF_FILE="${REFERENCE_DATA}/Homo_sapiens.GRCh38.100.gtf"
-FASTA_FILE="${REFERENCE_DATA}/Homo_sapiens.GRCh38.dna.primary_assembly.fa"
 
 STAR \
     --runThreadN 8 \
     --runMode genomeGenerate \
+    --readFilesCommand zcat \
     --sjdbOverhang 100 \
     --genomeFastaFiles ${FASTA_FILE} \
     --sjdbGTFfile ${GTF_FILE} \
@@ -50,9 +40,9 @@ STAR \
 # - We are assuming paired-end, gzip compressed (--readFilesCommand zcat)  FastQ
 #   files.
 
-# The compressed, paired-end, FastQ's after trimming (step 2).
+# THe compressed paired-end FastQ's that we are aligning.
-R1="${PROJECT_DIRECTORY}/step2/sample1_R1_paired.fastq.gz"
+R1="sample1_R1.fastq.gz"
-R2="${PROJECT_DIRECTORY}/step2/sample1_R2_paired.fastq.gz"
+R2="sample1_R2.fastq.gz"
 
 STAR \
     --runThreadN 8 \
@@ -60,4 +50,4 @@ STAR \
     --readFilesIn "${R1}" "${R2}" \
     --outSAMtype BAM SortedByCoordinate \
     --genomeDir ${GENOME_INDEX} \
-    --outFileNamePrefix "${ALIGNMENT_OUTPUT}/sample1_"
+    --outFileNamePrefix "${ALIGNMENT_OUTPUT}"

rnaseq/step4_multiqc/snippet.sh

@@ -1,5 +0,0 @@
-#!/bin/bash
-module load multiqc
-
-PROJECT_DIRECTORY="/groups/umcg-griac/tmp01/rawdata/$(whoami)/rnaseq"
-multiqc "${PROJECT_DIRECTORY}"

rnaseq/step5_count/snippet.sh

@@ -1,28 +0,0 @@
-#!/bin/bash
-PROJECT_DIRECTORY="/groups/umcg-griac/tmp01/rawdata/$(whoami)/rnaseq"
-COUNT_OUTPUT="${PROJECT_DIRECTORY}/step5"
-mkdir -p "${COUNT_OUTPUT}"
-
-# Storage location of annotation on Gearshift.
-REFERENCE_DATA="/groups/umcg-griac/prm03/rawdata/reference/genome"
-GTF_FILE="${REFERENCE_DATA}/Homo_sapiens.GRCh38.100.gtf"
-
-# Where our alignment file was stored.
-BAM="${PROJECT_DIRECTORY}/step3/alignment/sample1_Aligned.sortedByCoord.out.bam"
-
-# Compute counts using htseq-count.
-#
-# N.B.:
-# - If you are processing multiple files, consider using the `--nprocesses` flag
-#   to distribute computation of the files to different cores.
-# - The BAM file must be position sorted. If you used STAR with the
-#   `SortedByCoordinate` option you should be okay. If not, sort your BAM file
-#   using `samtools sort`.
-# - By default, strand aware library preparation is assumed. If not, specify the
-#   `--stranded` flag.
-htseq-count \
-    --order pos \
-    -f bam \
-    ${BAM} \
-    ${GTF_FILE} \
-    > ${COUNT_OUTPUT}/counts.txt

rnaseq/step6_overall_QC/01 - Principle Component Analysis.R

@@ -18,7 +18,7 @@ do.scale = FALSE
 
 
 # The analysis
-# We use prcomp to calculate the PCAs. Afterwards you should plot the results.
+# We ise prcomp to calculate the PCAs. Afterwards you should plot the results.
 norm.expr.data <- expression.data %>%
 	tibble::column_to_rownames("Gene")
 norm.expr.data <- norm.expr.data[rowSums(norm.expr.data) >= 10,] %>%
@@ -37,8 +37,8 @@ norm.expr.data.pcs <- norm.expr.data %>%
   )
 
 # Write summary of PCAs to files
-pcs.summary <- summary(norm.expr.data.pcs)
+pcs.summery <- summary(norm.expr.data.pcs)
-pcs.summary$importance %>%
+pcs.summery$importance %>%
 	t() %>%
 	as.data.frame() %>%
 	tibble::rownames_to_column("PC.name") %>%
@@ -46,7 +46,7 @@ pcs.summary$importance %>%
 		file.path(results.dir.pca, "importance.csv")
 	)
 
-pcs.summary$x %>%
+pcs.summery$x %>%
 	t() %>%
 	as.data.frame() %>%
 	tibble::rownames_to_column("ensembl.id") %>%
@@ -54,7 +54,7 @@ pcs.summary$x %>%
 		file.path(results.dir.pca, "values.csv")
 	)
 
-pcs.summary$rotation %>%
+pcs.summery$rotation %>%
 	t() %>%
 	as.data.frame() %>%
 	tibble::rownames_to_column("sample.id") %>%
@@ -63,15 +63,15 @@ pcs.summary$rotation %>%
 	)
 
 data.frame(
-		rownames = names(pcs.summary$center),
+		rownames = names(pcs.summery$center),
-		center = pcs.summary$center,
+		center = pcs.summery$center,
-		scale = pcs.summary$scale
+		scale = pcs.summery$scale
 	) %>%
 	readr::write_csv(
 		file.path(results.dir.pca, "rest.csv")
 	)
 
-# Not saved: pcs.summary$sdev,
+# Not saved: pcs.summery$sdev,
 
 
 # Next thing to do:

rnaseq/umi/07_remove_dups.pl

@@ -1,61 +0,0 @@
-#!/usr/bin/perl -w
-use strict;
-use Parallel::ForkManager;
-# this script creats one file with UMI unique reads and one with UMI duplicated reads
-# as input you need aligned sorted by coordinate bam file
-my @torun = ();
-foreach my $file ( <*Aligned.sortedByCoord.out.bam> ) {
-    push @torun, $file;
-}
-
-my $pm = Parallel::ForkManager->new( 12 );
-foreach my $file ( @torun ) {
-    my $sample = $file;
-    $sample =~ s/Aligned\.sortedByCoord\.out\.bam$//;
-    next if -s $sample.'_uniq.bam';
-    warn "Parsing $sample\n";
-    $pm->start and next;
-    my %seen = ();
-    my %duplicates = ();
-    my ( $uniqs, $dups ) = ( 0, 0 );
-    open F, 'samtools view -h '.$file.' |';
-    open F1, '| samtools view -bS - >'.$sample.'_uniq.bam';
-    open F2, '| samtools view -bS - >'.$sample.'_dups.bam';
-    while ( <F> ) {
-        if ( m/^\@/ ) {
-            print F1;
-            print F2;
-            next;
-        }
-        my ( $id, $flag, $chr, $pos, $mapq, $cigar, $chr2, $pos2, $tlen ) = split /\t/;
-        next if $flag & 256 or $flag & 512 or $flag & 1024; #skip if the read is not primary alignment/read fails platform/vendor quality checks/read is PCR or optical duplicate
-#        foreach ( 256, 512, 1024 ) { $flag-=$_ if $flag&$_ }
-        my ( $bc ) = $id =~ m/\:([GATCN\d]+)$/; #extract UMI barcode
-        my $uniq = join( ':', $chr, $pos, $flag, $tlen, $bc );
-        my $pos_ = $pos-1;
-        while ( $cigar =~ m/(\d+)([SHMDIN=])/g ) {
-            $pos_+=$1 if $2 eq 'M' or $2 eq '=' or $2 eq 'D' or $2 eq 'N'; # find position for minus strand
-        }
-        my $uniq2 = join( ':', $chr, $pos_, $flag, $tlen, $bc );
-        if ( exists($duplicates{$id}) or # already marked as duplicate
-            ( not($flag&16) and ++$seen{ $uniq } > 1 ) or # plus strand
-            ( $flag&16 and ++$seen{ $uniq2 } > 1 ) #minus strand
-            ) {
-            print F2;
-            $dups++;
-            $duplicates{$id} = 1;
-        }
-        else {
-            print F1;
-            $uniqs++;
-        }
-    }
-    close F;
-    close F1;
-    close F2;
-    warn "Unique: $uniqs (", 100*$uniqs/($uniqs+$dups), ")\n";
-    system( 'samtools', 'index', $sample.'_uniq.bam' );
-#    last;
-    $pm->finish;
-}
-$pm->wait_all_children;