Use BAM files for htseq-count.

Reviewed-on: #8

rnaseq/step2_trim/snippet.sh

@@ -33,7 +33,7 @@ java -jar $EBROOTTRIMMOMATIC/trimmomatic.jar PE \
   "${FASTQ_OUT}/sample1_R1_unpaired.fastq.gz" \
   "${FASTQ_OUT}/sample1_R2_paired.fastq.gz" \
   "${FASTQ_OUT}/sample1_R2_unpaired.fastq.gz" \
-  ILLUMINACLIP: TruSeq3-PE.fa:2:30:10 \
+  ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
   LEADING:3 \
   TRAILING:3 \
   SLIDINGWINDOW:4:25 \

rnaseq/step3_align/snippet.sh

@@ -2,33 +2,42 @@
 #
 # Align reads against reference genome.
 
+REFERENCE_DATA="/groups/umcg-griac/prm03/rawdata/reference/genome"
+
 PROJECT_DIRECTORY="/groups/umcg-griac/tmp01/rawdata/$(whoami)/rnaseq"
-# Store genome index in this location:.
-GENOME_INDEX="${PROJECT_DIRECTORY}/step3/genome_index"
-mkdir -p "${GENOME_INDEX}"
-# Store the generated `Aligned.sortedByCoord.out.bam` in this dir.
-ALIGNMENT_OUTPUT="${PROJECT_DIRECTORY}/step3/alignment/"
+# Store the generated `sample1_Aligned.sortedByCoord.out.bam` in this dir.
+ALIGNMENT_OUTPUT="${PROJECT_DIRECTORY}/step3/alignment"
 mkdir -p "${ALIGNMENT_OUTPUT}"
 
-# 1) Generate genome index.
+# 1) Generate genome index (optional).
+#
+# 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
 #   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.
-# - We're using gzip compressed reference data (--readFilesCommand zcat), i.e.,
-#   .gtf.gz and fa.gz. If not, you can remove the `zcat` flag.
+# - If you're using gzip compressed reference data, i.e., .gtf.gz and fa.gz,
+#   pass the `--readFilesCommand zcat` flag.
 
-# Storage location reference data (in this case on Gearshift).
-REFERENCE_DATA="/groups/umcg-griac/prm03/rawdata/reference/genome"
-GTF_FILE="${REFERENCE_DATA}/Homo_sapiens.GRCh38.100.gtf.gz"
-FASTA_FILE="${REFERENCE_DATA}/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz"
+# Store created genome index in this location:.
+GENOME_INDEX="${PROJECT_DIRECTORY}/step3/genome_index"
+mkdir -p "${GENOME_INDEX}"
+
+# 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} \
@@ -51,4 +60,4 @@ STAR \
     --readFilesIn "${R1}" "${R2}" \
     --outSAMtype BAM SortedByCoordinate \
     --genomeDir ${GENOME_INDEX} \
-    --outFileNamePrefix "${ALIGNMENT_OUTPUT}"
+    --outFileNamePrefix "${ALIGNMENT_OUTPUT}/sample1_"

rnaseq/step5_count/snippet.sh

@@ -0,0 +1,28 @@
+#!/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/umi/07_remove_dups.pl

@@ -2,6 +2,7 @@
 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;
@@ -33,12 +34,12 @@ foreach my $file ( @torun ) {
         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'; 
+            $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 )
+            ( $flag&16 and ++$seen{ $uniq2 } > 1 ) #minus strand
             ) {
             print F2;
             $dups++;