Dear all,
I want to quantify gene expression for my Trinity transcriptome assembly using RSEM. I am getting confusing results depending on how I specify quality scores of the Illumina paired end reads. They are Illumina 1.5 and have phred 64 quality scores.
The RSEM manual says that phred+33 is default:
--phred33-quals Input quality scores are encoded as Phred+33. (Default: on)
--phred64-quals Input quality scores are encoded as Phred+64 (default for GA Pipeline ver. >= 1.3). (Default: off)
a. When I specify the following commands I get the following alignment stats
rsem-calculate-expression --paired-end 1.out_paired.fastq 2.out_paired.fastq Trinity.fasta Trinity --phred64-quals --bowtie2 -p 32
19618532 reads; of these: 19618532 (100.00%) were paired; of these: 3453455 (17.60%) aligned concordantly 0 times 3684847 (18.78%) aligned concordantly exactly 1 time 12480230 (63.61%) aligned concordantly >1 times 82.40% overall alignment rate
bowtie2 -q --phred64 --sensitive --dpad 0 --gbar 99999999 --mp 1,1 --np 1 --score-min L,0,-0.1 -I 1 -X 1000 --no-mixed --no-discordant -p 10 -k 200 -x Trinity. fasta -1 1.out_paired.fastq -2 2.out_paired.fastq | samtools view -S -b -o temp/ bam -
b. However, I get exactly the same stats when I don’t specify phred 64. How can this be if RSEM defaults to phred33?
rsem-calculate-expression --paired-end 1.out_paired.fastq 2.out_paired.fastq Trinity. fasta Trinity --bowtie2 -p 32
19618532 reads; of these: 19618532 (100.00%) were paired; of these: 3453456 (17.60%) aligned concordantly 0 times 3684845 (18.78%) aligned concordantly exactly 1 time 12480231 (63.61%) aligned concordantly >1 times 82.40% overall alignment rate
I checked in the RSEM standard output and it is specifying phred33. So phred33 is default but how can it give the same alignment?
bowtie2 -q --phred33 --sensitive --dpad 0 --gbar 99999999 --mp 1,1 --np 1 --score-min L,0,-0.1 -I 1 -X 1000 --no-mixed --no-discordant -p 10 -k 200 -x Trinity. fasta -1 1.out_paired.fastq -2 2.out_paired.fastq | samtools view -S -b -o temp/ bam -
c. When I directly specify phred33 I also get the same alignment stats
rsem-calculate-expression --paired-end 1.out_paired.fastq 2.out_paired.fastq Trinity.fasta Trinity --phred33-quals --bowtie2 -p 32
19618532 reads; of these: 19618532 (100.00%) were paired; of these: 3453456 (17.60%) aligned concordantly 0 times 3684845 (18.78%) aligned concordantly exactly 1 time 12480231 (63.61%) aligned concordantly >1 times 82.40% overall alignment rate
Any help would be much appreciated.
A
