Paired end inconsistency in tophat SAM

[dariober]

Hi there,

Apologies if I missing something obvious here...
In the SAM output of tophat (accepted_hits.sam) I found read pairs where the two reads mapped to different chromosomes (3rd column). However the column MRNM (mate reference sequence, 7th column) has '=' meaning both reads mapped to the same reference.

Here's an example:

EBRI093151_0001:8:100:268:457#NNNNNN	145	10	65535404	255	35M	=	85418059	0	ATTTGCACTATTACACTTAAATTGTTATCCTTTTT	_aaa`b^baaaa_b_baabaaababababbbabba	NM:i:2
EBRI093151_0001:8:100:268:457#NNNNNN	97	5	85418059	255	35M	=	65535404	0	ATTTTTACATCAGTCCCTTTAACACAAATCCATAT	aabaabaabaaaa^a_a`a`^``a_a_aa_^]`aa	NM:i:0

I discovered this while running HTSeq to count reads mapping to GTF features. I get the warning:

Warning: Incorrect 'proper_pair' flag value for read pair EBRI093151_0001:8:100:269:1653#NNNNNN

For info, this is the tophat command I used:

tophat \
  --output-dir /exports/.../20100122_RNAseq_CTRL/ \
  --mate-inner-dist 130 \
  --mate-std-dev 30 \
  --solexa1.3-quals \
  --GFF /exports/.../Sus_scrofa.Sscrofa9.56.gff3 \
  /exports/.../bowtie/current/indexes/Sscrofa9.56_plus_Human_ribosomal_DNA_complete_repeating_unit \
  /exports/.../RNAseq_CTRL_1_35bp.fastq \
  /exports/.../RNAseq_CTRL_2_35bp.fastq

Is this a bug in tophat (tophat-1.0.13)? Does anyone know whether such pairs are reliable?

Many thanks

Dario

[Daehwan]

Hi there,

Apologies if I missing something obvious here...
In the SAM output of tophat (accepted_hits.sam) I found read pairs where the two reads mapped to different chromosomes (3rd column). However the column MRNM (mate reference sequence, 7th column) has '=' meaning both reads mapped to the same reference.

Here's an example:

EBRI093151_0001:8:100:268:457#NNNNNN	145	10	65535404	255	35M	=	85418059	0	ATTTGCACTATTACACTTAAATTGTTATCCTTTTT	_aaa`b^baaaa_b_baabaaababababbbabba	NM:i:2
EBRI093151_0001:8:100:268:457#NNNNNN	97	5	85418059	255	35M	=	65535404	0	ATTTTTACATCAGTCCCTTTAACACAAATCCATAT	aabaabaabaaaa^a_a`a`^``a_a_aa_^]`aa	NM:i:0

I discovered this while running HTSeq to count reads mapping to GTF features. I get the warning:

Warning: Incorrect 'proper_pair' flag value for read pair EBRI093151_0001:8:100:269:1653#NNNNNN

For info, this is the tophat command I used:

tophat \
  --output-dir /exports/.../20100122_RNAseq_CTRL/ \
  --mate-inner-dist 130 \
  --mate-std-dev 30 \
  --solexa1.3-quals \
  --GFF /exports/.../Sus_scrofa.Sscrofa9.56.gff3 \
  /exports/.../bowtie/current/indexes/Sscrofa9.56_plus_Human_ribosomal_DNA_complete_repeating_unit \
  /exports/.../RNAseq_CTRL_1_35bp.fastq \
  /exports/.../RNAseq_CTRL_2_35bp.fastq

Is this a bug in tophat (tophat-1.0.13)? Does anyone know whether such pairs are reliable?

Many thanks

Dario

The current version (1.1.2) checks the distance between the paired reads. If the distance is longer than #mate-inner-dist + mate-std-dist#, it doesn't output such a pair. So, I think it's probably a bug in the version (1.0.13) you use. Can you please run it again with 1.1.2 or we'll release 1.1.3 soon with some fix in strand-specific RNA-Seq and some additions.

[dariober]

The current version (1.1.2) checks the distance between the paired reads.

Thanks Daehwan, I will switch to 1.1.2. By the way, in some post, which I can't find again, I think I read that tophat 1.1.x doesn't require to specify the options --mate-inner-dist and --mate-std-dist since it will infer them by its own. Is this correct?

If the distance is longer than #mate-inner-dist + mate-std-dist#, it doesn't output such a pair.

Isn't a distance of "mate-inner-dist + mate-std-dist" too conservative to reject a 'proper pair'? I would say that something like mate-inner-dist + 2(or 3)stdev is still acceptable (just a thought...).

So, I think it's probably a bug in the version (1.0.13) you use. Can you please run it again with 1.1.2 or we'll release 1.1.3 soon with some fix in strand-specific RNA-Seq and some additions.

Looking forward to new updates then!

All the best
Dario