Hello! I analyzing paired-end RNA-seq reads made from a strand-specific dUTP Illumina library construction kit (Read2 corresponding to the original transcript orientation). These reads are from one organism that has a reference genome, and I am hoping to perform differential expression analyses with HT-seq. I have aligned my strand-specific reads to the genome with Bowtie2 (using the --fr --nofw flags), and now I need to pull the appropriate alignments out of the SAM output file.
I've read Istvan Albert's great post on using SAMtools to split strand-specific data, and looked at how SAM flags are visualized, but I am still lost. My original thought was to output SAM flag -3 (read paired, read mapped in proper pair), but does this not take in to consideration strand-specificity? What are the appropriate flags to pull out proper mate-pairs of reads that also correspond to the strand that it originated from? Is this extremely important for downstream HT-seq counts to get differential expression?
Thank you!
•
written
by
Etstevens •
2 answers
Hi a few things. There's a bit of confusion so I;m going to clear things up. I hope you don't mind the corrections.
differential expression analyses with HT-seq
HTSeq count script does not perform DE analysis. That's for tools like edgeR or DESeq2.
Bowtie2 (using the --fr --nofw flags)
Why are you excluding unpaired reads? IMO this is not needed for DE analysis. Also, just to confirm your organism is for example a bacterium that does not undergo mRNA splicing? If it does, you need to use a splice aware aligner such as STAR or HISAT2. Bowtie2 is not splice aware.
Is this extremely important for downstream HT-seq counts to get differential expression?
From my understanding no and... yes, it's complicated. You don't have to worry about filtering out non-stranded non-mate read pairs. HTSeq from my understanding uses a (relatively) simple counting method described here. Depending on what you select as the count method, it will change how unaligned reads are counted.
BTW, a very important flag is --stranded option. You should specify --stranded yes or --stranded reverse. reverse is for some illumina library prep that produces reads in reverse orientation (so forward = reverse, etc). Unsure what which your library is. Also you need to sort your sam files by read name before counting.
The reason I say no and yes is because you need a good reason for filtering out non-paired or non stranded reads. Anything performed upstream in the alignment/counting can have major ramification downstream in the DE analysis. From my reading the best approaches use very little filtering/trimming of the reads upstream. Whether a paired/non-paired should be counted/filtered is a question I don't know the answer to. But IMO you should count non-paired reads.
I hope that answers your questions.
I agree with @Amar - there is no need for requiring the paired end information but standness is important. But why not use some of the more modern pseudo aligners for mappung/quantification? You can read about such considerations here - they also give better DE results as I answered in this previous answer.
Log in to answer this question.
-
Same strand, same QNAME different FLAG?
written by compuTE
-
Which column to use from STAR output for Differential Gene Expression Analysis?
written by caggtaagtat
-
Choosing sam flags to extract all reads mapped to 5' -> 3' OR 3' -> 5' of reference sequence?
written by johnnytam100
-
Do you use anti-sens reads for differential gene expression analysis with Stranded Paired-End Libra…
written by ZheFrench
-
Single-end strand specific Rna-Seq: how to separate alignments by strand
written by termanini •
-
Samtools mpileup vcf file
written by AW
-
strandedness - RNASeq library preparation
written by Sashi Kiran Challa
-
Paired-end stranded RNA-seq...still confused about TopHat
written by cadeans •
-
SAM flags for strand-specific protocols
written by Pfs
-
Ht-Seq Read Count And Strand-Specificity
written by thecuriousbiologist
