search

DaehwanKimLab / tophat

Spliced read mapper for RNA-Seq
http://ccb.jhu.edu/software/tophat
Boost Software License 1.0
90 stars 46 forks source link

Tophat v2.1.1 reports identical alignments and does not mark either as secondary. #36

Open MatthewRalston opened 8 years ago

MatthewRalston commented 8 years ago

Tophat 2.1.1 reports identical alignments, related to a similar issue with HISAT2. The SAM records are identical, and cause troubles with downstream tools. Is there a way to limit these duplicate alignments? Or should we wait until this is fixed

Sofware versions: samtools: 1.3.1 Picard 2.2.2 Tophat 2.1.1

 $TOPHAT2_DIR/tophat2 -p $THREADS \
            --no-coverage-search \
            --rg-id 1 --rg-sample $SAMPLE --rg-library $SAMPLE --rg-platform illumina \
            $REFERENCE_DIR/$REFERENCE/bowtie2_index/$REFERENCE \
            $SAMPLE.uncontaminated.fastq.1.gz $SAMPLE.uncontaminated.fastq.2.gz

$SAMTOOLS merge - tophat_out/accepted_hits.bam tophat_out/unmapped.bam | $SAMTOOLS sort - mark_dup                                    
 ${JAVA} -Xmx4G -jar $PICARD_JAR MarkDuplicates \
             INPUT=mark_dup.bam \
             OUTPUT=/dev/null \
             METRICS_FILE=${SAMPLE}.duplicate_stats.txt \
             VALIDATION_STRINGENCY=SILENT;

Picard Exception:

Exception in thread "main" htsjdk.samtools.SAMException: Value was put into PairInfoMap more than once.  96: null:NB501257:32:HWCMVBGXX:1:21209:24100:6206
        at htsjdk.samtools.CoordinateSortedPairInfoMap.ensureSequenceLoaded(CoordinateSortedPairInfoMap.java:133)
        at htsjdk.samtools.CoordinateSortedPairInfoMap.remove(CoordinateSortedPairInfoMap.java:86)
        at picard.sam.markduplicates.util.DiskBasedReadEndsForMarkDuplicatesMap.remove(DiskBasedReadEndsForMarkDuplicatesMap.java:61)
        at picard.sam.markduplicates.MarkDuplicates.buildSortedReadEndLists(MarkDuplicates.java:442)
        at picard.sam.markduplicates.MarkDuplicates.doWork(MarkDuplicates.java:193)
        at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:209)
        at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:95)
        at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:105)

Locating the offending reads:

samtools view mark_dup.bam | grep NB501257:32:HWCMVBGXX:1:21209:24100:6206
NB501257:32:HWCMVBGXX:1:21209:24100:6206        177     chr1    113980481       50      75M     chr12   95825708        0       TCTTCTTCTTCTTCTTCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTCCTTCTCCT  EEEEEEEEEEAAEEAEEAEEE/EEEEE/EEEEE/EEEEE/EEEEEAEEEEEEEAEEEEEEEEEEEEEE6EAAAA6 AS:i:-5  XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:5C69       YT:Z:UU NH:i:1
NB501257:32:HWCMVBGXX:1:21209:24100:6206        177     chr1    113980481       50      75M     chr12   95825708        0       TCTTCTTCTTCTTCTTCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTCCTTCTCCT  EEEEEEEEEEAAEEAEEAEEE/EEEEE/EEEEE/EEEEE/EEEEEAEEEEEEEAEEEEEEEEEEEEEE6EAAAA6 AS:i:-5  XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:5C69       YT:Z:UU NH:i:1
NB501257:32:HWCMVBGXX:1:21209:24100:6206        113     chr12   95825708        50      76M     chr1    113980481       0       AGAAGTAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAA AE/EE/EE/E/EEA/EE<AE/EA/EAEEAEEEEEEEEEEE<AEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAAAS:i:0   XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:76 YT:Z:UU NH:i:1
NB501257:32:HWCMVBGXX:1:21209:24100:6206        113     chr12   95825708        50      76M     chr1    113980481       0       AGAAGTAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAAGAA AE/EE/EE/E/EEA/EE<AE/EA/EAEEAEEEEEEEEEEE<AEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAAAS:i:0   XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:76 YT:Z:UU NH:i:1