Entering edit mode
13.0 years ago
Serena
•
0
Hi guys, I am currently working on paired-end RNA-seq data. I am running samtools using the following script for single ends, but I am not sure it is correct for paired-end, specially the part I use to get uniq.bam file.
samtools view -bS s1sequence.sam.gz > s1sequence.bam
samtools view -bq 1 s1sequence.bam > s1uniq.bam
samtools sort s1uniq.bam s1uniq.sorted
samtools index s1uniq.sorted.bam
samtools mpileup -vcf wg.fa s1uniq.sorted.bam > s1.pileupraw
Do you know if it is correct for paired-ends? I am a beginner so any advise is more than welcome! Thank you very much in advance!
Serena
What do you want uniq.bam to contain?
looks fine, you can combine the first 2 steps with: samtools view -bSq 1 sequence.sam.gz > s1unique.bam
Sort of a genreal question...is mapping quality independant for both reads? Or would it take into account if the mate mapped uniquely? Because I would think that RNASeq would have a lot of pairs where one end fell in a repeated domain, but the other end was unique.
-v and -c are not mpileup options, those are old pileup options
Are you sure you want to run mpileup on the data - do you want SNP's/indels or per-base sequencing depth for an RNA-seq dataset?
Hi! thank you very much for your comments and help! @Sean, As for the uniq.bam, it should contain unique mapped reads @swbarnes2 Yes, I would like to take into account just the mate mapped uniquely...but I don't know how to do that! @Jeremy Thanks a lot for telling me about the old pileup options! @Chris Penkett, Yes,I would like to pileup SNPs/indels. Again thank everybody in advance for your help! Serena
errata corrige. uniq bam should contain mapped reads where one end map uniquely and the other ambiguously. Is correct to use the following script to get that? samtools view -bq 1 s1sequence.bam > s1uniq.bam thanks for your help in advance!