如何對多個轉錄組測序數據找變異呢

• 第一次對參考基因組建索引

• 然後進行第一次序列比對

• 之後根據第一次比對得到的所有剪切位點，重新對參考基因組建立索引

• 再進行 STAR 二次序列比對

由於後面要用 GATK 進行 call 變異，還需要對比對結果 SAM 文件進行一些處理，這些都可以用 picard 來做，包括 使用AddOrReplaceReadGroups命令SAM頭文件添加 @RG 標籤，SAM 文件排序並轉 BAM 格式，然後使用MarkDuplicates命令標記 duplicate，由於 STAR 軟體使用的 MAPQ 標準與 GATK 不同，而且比對時會有 reads 的片段落到內含子區間，需要進行一步 MAPQ 同步和 reads 剪切，使用 GATK 專為 RNA-seq 應用開發的工具 進行操作，它會將落在內含子區間的 reads 片段直接切除，並對 MAPQ 進行調整。

問題描述：

GATK Best Practices for RNA-seq suggests using the STAR 2-pass alignment steps to do the mapping. It sounds very good. But in operation, maybe there are some issues. As shown in the methods, for every sample, the genome index should be generated twice, which is time-consuming for big genome such as human. Moreover, it seems that we can not do alignment for two samples at same time with the same genome path, right? Should we copy a genome.fa for every sample in an individual place? Seems crazy for a big project.

解答：

其實可以首先把fastq批量比對後得到多個.tab後綴文件，直接cat到一起作為第二次構建index的參考即可。

比對後得到的bam文件批量走gatk的找變異流程，代碼如下：

喜歡這篇文章嗎？立刻分享出去讓更多人知道吧！