零基础入门转录组上游分析——第四章（序列比对）

零基础入门转录组上游分析——第四章序列比对目录零基础入门转录组上游分析——第四章序列比对1. 之前章节结果的查看1. 构建参考基因组索引2. 序列比对3. 压缩和排序XXX.sam文件4. 构建bam文件的索引可选5. 进阶我这里使用的虚拟机是vmwarewokstation版本16.0.0 linux系统是ubantu64位版本20.04.3 上一章我们做了质控分析并对数据进行了质量过滤这一章序列比对将会用到质量过滤后的数据02_clean_data以及第二章准备的参考基因组文件。 本实验选用的是模式生物——C57BL/6J小鼠。 实验分组药物处理组对照组每组6只鼠。 软件hisat2, samtools1. 之前章节结果的查看在前两章中分别做了原始数据的准备和质控以及质量过滤。首先来看一下现有的文件和之前的输出结果如下图所示质控的结果我给删了有点占地方00_raw_data文件夹中是我们的原始数据和样本信息表但是有了clean_data这些就用不到了。01_ref文件夹中有参考基因组XXX.fa.gz和注释文件XXX.gtf.gz02_clean_data文件夹中有上一章做完质量过滤的数据又叫做clean data。在02_clean_data文件夹中有两个XXX.fq.gz文件由于我这里是双端测序所以 有A2_1_val_1, A2_2_val_2这两个。1. 构建参考基因组索引有了clean_data和准备好的参考基因组我们接下来第一步就是构建参考基因组索引要用到的软件是hisat2如果没有安装过这个软件的同学可以参考软件的安装。1先将01_ref文件夹中的参考基因组压缩包(XXX.fa.gz)和注释文件压缩包(XXX.gtf.gz)解压gunzip Mus_musculus.GRCm39.dna.primary_assembly.fa.gz解压参考基因组gunzip Mus_musculus.GRCm39.108.gtf.gz解压参考基因组注释文件解压后的结果如下图所示xxx.fa文件就是我们要用到参考基因组2接下来是构建参考基因组索引文件先通过pwd指令查看当前文件夹的路径接下来用hisat2-build 构建参考基因组的索引我这里用的是绝对路径hisat2-build /home/daizhuer/Desktop/01_ref/Mus_musculus.GRCm39.dna.primary_assembly.fa /home/daizhuer/Desktop/01_ref/genome 1hisat2-build.log 21hisat2-build指令不难hisat2-build后跟的是参考基因组的路径再后面是输出结果的路径我这里是输出到当前文件夹下命名为genome…再后面跟的是hisat2-build结果输出报告hisat2-build.log 文件构建参考基因组索引这步会很慢需要耐心等待结果如下图所示多了一些白色的genome…的文件这些就是构建好的参考基因组索引文件注意备份2. 序列比对构建好索引文件后接下来就要将过滤后的原始数据比对到参考基因组上首先我们先了解一下常用的参数-new-summary这会让hisat2软件比对后输出的报告会更加好看-p是线程数我的虚拟机是内存8G分配了4个处理器256G硬盘空间线程数跟处理器的数量有关这个线程数设的越大运行速度会越快根据你的虚拟机状况选用合适的线程数。-x这里输入构建好的参考基因组索引的路径-1输入过滤后的原始数据1的路径-2输入过滤后的原始数据2的路径-S输出的XXX.sam文件-rna-strandness这个参数是链特异性文库需要的就是说你要问测序公司你的原始数据是普通文库数据还是链特异性文库。如果是普通文库这里就不加-rna-strandness这个参数hisat2软件会比对两次正向一次反向互补的链再比对一次如果是链特异性文库就要加上-rna-strandness参数RF参数代表的是双端测序如果是单端测序这里就是R首先在桌面路径下通过mkdir 03_hisat2_result创建一个03_hisat2_result文件夹(用来存放比对结果)接下来看代码(我这里用的全是绝对路径绝对路径好处就是你在任何路径下运行代码都没问题但是缺点就是特别长会看花眼)hisat2 --new-summary -p 2 -x /home/daizhuer/Desktop/01_ref/genome -1 /home/daizhuer/Desktop/02_clean_data/A2_1_val_1.fq.gz -2 /home/daizhuer/Desktop/02_clean_data/A2_2_val_2.fq.gz -S /home/daizhuer/Desktop/03_hisat2_result/A2.sam --rna-strandness RF 1/home/daizhuer/Desktop/03_hisat2_result/A2.log 21注意hisat2比对生成的XXX.sam文件会非常大大概20G左右所以一定要留出足够的空间运行结果如下我们可以看到生成了两个文件一个是20G的A2.sam文件另一个是A2.log文件这里存放比对结果的分析报告。通过cat A2.log指令可以查看比对结果分析报告如下图所示结果显示97.49%的reads都比对到参考基因组了说明比对结果非常好3. 压缩和排序XXX.sam文件根据前面的图我们可以看到sam文件非常大这么大的数据不利于保存和分析现需要对其进一步压缩和排序。这里我们需要用到samtools软件如果没安装的小伙伴可以通过conda install samtools -y指令安装即可samtools软件可以轻松对sam文件进行压缩和排序输出结果为XXX.bam文件指令如下sort 参数表示进行排序-o 是输出文件的路径samtools sort -o /home/daizhuer/Desktop/03_hisat2_result/A2.bam /home/daizhuer/Desktop/03_hisat2_result/A2.sam代码大致意思就是输入03_hisat2_result文件夹中的A2.sam文件通过samtools软件对其压缩和排序输出结果到03_hisat2_result文件夹下命名为A2.bam。结果如下图所示压缩和排序后生成了一个1.7G的A2.bam文件相比于sam文件小了10倍有了bam文件sam文件就没啥用了可以删除节省空间如果电脑或移动硬盘比较大的同学可以选择两个都备份。注意XXX.log文件不要删这可以帮助后期查看比对率4. 构建bam文件的索引可选有了bam文件接下来可以对bam文件构建索引构建索引是为了能够在IGV软件中观察比对结果这步不是必须的可以根据需求选择性观看。构建bam文件的索引非常简单通过samtools index A2.bam指令就能完成。结果如下生成了一个2.2Mb的XXX.bai文件这就是索引文件通常是在IGV软件中使用可视化比对结果。5. 进阶在这一部分主要介绍一下如何批量比对1准备样本信息表切换到02_clean_data文件夹路径下输入指令ls *fq.gz sample_ifo会在当前路径下生成一个sample_ifo文件接下来通过vim sample_ifo指令编辑文件如下图所示是刚开的样子点一下键盘上的i键进入编辑模式修改成下图所示的样子中间以空格相隔一共三列前两列就是clean_data的名称第三列是你想要的样本名可以随便取自己记得是啥就行记住这三列等会会用到编辑好后摁下Esc输入:wq保存并退出。2脚本编辑有了样本信息表之后我们就可以开始写脚本了先切回到03_hisat2_result文件夹目录下输入如下指令awk{print hisat2 --new-summary -p 2 -x /home/daizhuer/Desktop/01_ref/genome -1 /home/daizhuer/Desktop/02_clean_data/$1 -2 /home/daizhuer/Desktop/02_clean_data/$2 -S /home/daizhuer/Desktop/03_hisat2_result/$3.sam --rna-strandness RF 1/home/daizhuer/Desktop/03_hisat2_result/$3.log 21}/home/daizhuer/Desktop/02_clean_data/sample_iforun_hisat2.sh输出为一个名为run_hisat2.sh的脚本文件代码较长你们可以粘贴下来慢慢理解其实和序列比对那里的代码很像只不过有一些需要用到名称的地方用刚才准备的样本信息表中不同列的名称给代替了。例如hisat2 --new-summary -p 2 -x /home/daizhuer/Desktop/01_ref/genome -1 /home/daizhuer/Desktop/02_clean_data/$1这个$1就是样本信息表中的第一列也就是刚才样本信息表中的A2_1_val_1.fq.gz后面的$2…依次类推如果你的样本信息表有多行那么生成的这个XXX.sh脚本每一行都有这样相似代码只不过分析的样本名称不同。脚本生成之后就可以通过bash XXX.sh命令运行了不过这里一旦运行就要保证有足够的空间因为一次比对就会生成20g的文件需要保证足够的空间。压缩和排序的自动化也类似用awk指令即可拼接出一个简单脚本。awk ‘{print “…”}’ 样本信息表 XXX.sh注意如果这里理解的比较困难可以查看第二章质控进阶那部分那里awk代码较短更容易理解。注意samtools软件安装后如果没法使用可以删除conda环境后重新安装samtools软件代码如下2023年7.17更新conda install -c bioconda samtools1.9 -y参考链接解决samtools: error while loading shared libraries: libcrypto.so.1.0.0/libncurses.so.5的问题结语以上就是零基础入门转录组上游分析——第四章序列比对的所有过程如果有什么需要补充或不懂的地方大家可以私聊我或者在下方评论。如果觉得本教程对你有所帮助希望广大学习者能够点赞收藏加关注关于我们我们的团队是领航生信如果大家想要系统学习常规SCI生信套路和流程或者了解更多生信相关知识可以在下方公众号链接找到我们~~~祝大家能够开心学习轻松学习在学习的路上少一些坎坷~~~目录部分跳转链接零基础入门生信转录组数据分析——导读