输入读数为 5500 万，但只有 100 万用于对齐

Question

[U]我使用 tophat (v2.1.0) 运行此代码，以使用来自基因组 (Homo_sapiens_UCSC_hg19) 的 bowtie2 基因组.bt2 索引从我的 RNA-seq fastq 文件中对齐读取 (bowtie2 (v2.2.6.0)) ( [/U]：

tophat2 -p 8 -G /home/ajsn6c/Desktop/Kumar_RNA-seq/Homo_sapiens_UCSC_hg19 /Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/hg19.gtf /home/ajsn6c/Desktop/Kumar_RNA-seq/Homo_sapiens_UCSC_hg19/Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/genome HPDE_S11_L002_R1_001.fastq

[U]我的 fastq 文件大约 13 GB。但是，对齐后我接受的命中文件只有 50 MB。[/U]

[U]这里的对齐输出说我有大约 5500 万个保持读取：[/U]

[2018-02-21 13:58:33] 开始 TopHat 运行 (v2.1.0)

[2018-02-21 13:58:33]     Checking for Bowtie
      Bowtie version:    2.2.6.0
[2018-02-21 13:58:33] Checking for Bowtie index files (genome)..
[2018-02-21 13:58:33] Checking for reference FASTA file
[2018-02-21 13:58:33] Generating SAM header for /home/ajsn6c/Desktop /Kumar_RNA-seq/Homo_sapiens_UCSC_hg19/Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/genome
[2018-02-21 13:58:35] Reading known junctions from GTF file
[2018-02-21 13:58:39] Preparing reads
 left reads: min. length=12, max. length=101, 55970267 kept reads (45104 discarded)
Warning: short reads (<20bp) will make TopHat quite slow and take large amount of memory because they are likely to be mapped in too many places
[2018-02-21 14:17:45] Building transcriptome data files Panc1/tmp/genes
[2018-02-21 14:17:59] Building Bowtie index from genes.fa
[2018-02-21 14:32:14] Mapping left_kept_reads to transcriptome genes with Bowtie2 
[2018-02-21 15:38:44] Resuming TopHat pipeline with unmapped reads
[2018-02-21 15:38:44] Mapping left_kept_reads.m2g_um to genome genome with Bowtie2 
[2018-02-21 16:17:07] Mapping left_kept_reads.m2g_um_seg1 to genome genome with Bowtie2 (1/4)
[2018-02-21 16:18:13] Mapping left_kept_reads.m2g_um_seg2 to genome genome with Bowtie2 (2/4)
[2018-02-21 16:19:32] Mapping left_kept_reads.m2g_um_seg3 to genome genome with Bowtie2 (3/4)
[2018-02-21 16:20:46] Mapping left_kept_reads.m2g_um_seg4 to genome genome with Bowtie2 (4/4)
[2018-02-21 16:21:59] Searching for junctions via segment mapping
[2018-02-21 16:25:24] Retrieving sequences for splices
[2018-02-21 16:27:18] Indexing splices
Building a SMALL index
[2018-02-21 16:27:37] Mapping left_kept_reads.m2g_um_seg1 to genome segment_juncs with Bowtie2 (1/4)
[2018-02-21 16:27:50] Mapping left_kept_reads.m2g_um_seg2 to genome segment_juncs with Bowtie2 (2/4)
[2018-02-21 16:28:03] Mapping left_kept_reads.m2g_um_seg3 to genome segment_juncs with Bowtie2 (3/4)
[2018-02-21 16:28:17] Mapping left_kept_reads.m2g_um_seg4 to genome segment_juncs with Bowtie2 (4/4)
[2018-02-21 16:28:31] Joining segment hits
[2018-02-21 16:31:02] Reporting output tracks

[2018-02-22 19:21:42] A summary of the alignment counts can be found in ./tophat_out/align_summary.txt
[2018-02-22 19:21:42] Run complete: 02:08:37 elapse

[U]这是来自 align_summary 文件的对齐摘要[/U]：

reads:
      Input     :    926337
       Mapped   :    898584 (97.0% of input)
        of these:     14621 ( 1.6%) have multiple alignments (14 have >20)

97.0% 的总体读取映射率。

为什么输入只有 900K，而它保持了 5500 万次读取？读取的质量也具有出色的 phred 分数。任何想法将不胜感激！

谢谢 亚历克斯

Answer 1

您的日志文件中的这些条目很奇怪：

[2018-02-21 14:17:45] 构建转录组数据文件Panc1/tmp/genes

[2018-02-21 14:17:59] 从genes.fa构建Bowtie索引

这是您的 tophat2 命令（我已重新构建命令以提高可读性）

./tophat2 \
    -p 8 \
    -G /home/ajsn6c/Desktop/Kumar_RNA-seq/Homo_sapiens_UCSC_hg19 /Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/hg19.gtf \
    /home/ajsn6c/Desktop/Kumar_RNA-seq/Homo_sapiens_UCSC_hg19/Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/genome \
    HPDE_S11_L002_R1_001.fastq

1. 似乎有一些错误的空格（例如[...]Homo_sapiens_UCSC_hg19 /Homo_sapiens[...]；不确定这是否是问题所在。
2. 根据您的命令，转录组应基于文件[...]/UCSC/hg19/Sequence/Bowtie2Index/hg19.gtf 构建；我不知道Panc1/tmp/genes 来自哪里，但显然这个文件是用来构建参考转录组的，而不是[...]/hg19.gtf。