如何将测序reads比对至两条短参考序列并统计比对占比？

一、推荐比对工具 #

优先用Bowtie2，其次是BWA-MEM。这俩都是Illumina短reads比对的主流工具，针对你的150bp双端reads和短扩增子参考序列，速度快且比对准确性高，能自动选择每条reads的最佳比对结果，正好满足你区分两条参考序列的需求。

二、具体操作步骤 #

1. 整理参考序列 #

把两条参考序列合并成一个fasta文件（比如命名为virus_amplicons.fasta），给每条序列加上唯一ID，方便后续统计：

>virus1
aaaaagtataaatataggaccaggcagagcattttatacaacaggagaaataataggagatataagacaagcacattgtaaccttagtagagcaaaatggaatgacactttaaataagatagttataaaattaagagaacaatttgggaataaaacaatagtctttaagcact
>virus2
aaaaagtatccgtatccagaggggaccagggagagcatttgttacaataggaaaaataggaaatatgagacaagcacattgtaacattagtagagcaaaatggaatgccactttaaaacagatagctagcaaattaagagaacaatttggaaataataaaacaataatctttaagcaat

2. 构建索引 #

用Bowtie2构建索引（如果选BWA就用bwa index命令）：

bowtie2-build virus_amplicons.fasta virus_amplicons

3. 双端reads比对 #

对单个样本执行比对，输出SAM格式文件（加--no-unal只保留比对上的reads，节省空间）：

bowtie2 -x virus_amplicons -1 your_sample_R1.fastq.gz -2 your_sample_R2.fastq.gz --no-unal -S sample_aln.sam

4. 统计比对到两条参考序列的reads数 #

用samtools筛选出成功比对的reads，再分别统计对应两条病毒序列的数量：

# 统计比对到virus1的成对reads数（-f 2表示双端都比对成功）
samtools view -f 2 sample_aln.sam | grep -c "virus1"

# 统计比对到virus2的成对reads数
samtools view -f 2 sample_aln.sam | grep -c "virus2"

如果不想只统计成对比对的，去掉-f 2即可，但建议用成对统计更准确，避免单端假阳性。

5. 计算占比 #

总有效比对reads对数 = virus1数 + virus2数
virus1占比 = (virus1数 / 总有效数) * 100%
virus2占比 = (virus2数 / 总有效数) * 100%

三、额外建议 #

• 比对前先做质控：用fastp去除接头、截掉低质量末端，避免垃圾reads干扰比对结果。
• 多样本批量处理：写个shell循环脚本，自动遍历所有样本的fastq文件，完成比对和统计，省得手动重复操作。
• 若遇到少量reads同时匹配两条参考序列，工具会自动分配到比对质量更高的那条，不用额外处理，直接统计就行。

内容的提问来源于stack exchange，提问作者Sara Nicholson