[Linux|生信]project4_03:序列比对

背景

前面介绍了如何从sra文件得到fastq文件，并利用fastqc、multiqc生成质控报告以及使用trim_galore去除低质量碱基

[Linux|生信]project4_01:批量下载sra文件并转化为fastq文件

[Linux|生信]project4_02:质控过滤

今天介绍一下如何自建序列索引以及序列比对

Just Do It!

输入文件准备

1. 基因组注释文件（*.gff或*.gft）
2. 参考基因组文件（*.fa或*.fa.gz）
3. 前述经过质控的序列文件( .fa或.fa.gz） 注：本文后面会介绍如何下载基因组注释文件以及参考基因组文件

自建索引（hisat2）

# 自建索引（hisat2）
## 其实hisat2-build在运行的时候也会自己寻找exons和splice_sites，但是先做的目的是为了提高运行效率
extract_exons.py ../GRCh38_latest_genomic.gff > hg38.exons.gtf
extract_splice_sites.py ../GRCh38_latest_genomic.gff > hg38.splice_sites.gtf
# 建立index， 必须选项是基因组所在文件路径和输出的前缀
hisat2-build -p 16 --ss hg38.splice_sites.gtf --exon hg38.exons.gtf ../GRCh38_latest_genomic.fna hg38

生成的索引文件

参数说明

• extract_exons.py与extract_splice_sites.py此两个Python脚本为hisat2软件内置的脚本，成功安装即可使用
• 参考基因组与基因组注释文件均位于当前目录的上一级目录，操作时注意修改相应文件路径
• hisat2-build用法介绍

• -p: 指定计算的线程数，可加速运算
• --ss 和 --exon 后接文件为前述生成 .fna为参考基因组文件，正常为.fa。我下载的为*.fna格式，打开看了一下，里面存在大量的不确定碱基被表注为“N”，猜想*.fna文件为包含不确定碱基“N”的参考基因组文件。不影响后续操作。 最后的hg38为自定义的索引前缀，运行完发现生成了8个索引文件

比对

质控后得到的三对双端序列文件

• 首先将质控后的双端序列文件拷贝至一个空文件夹，以确保其他*.fq.gz文件不会对结果造成干扰
• 下述所有操作均在该文件夹中进行

处理思路

hisat2 -p 2 -x $HOME/Project4/source_data/cleandata/refANDgtf/index_by_hisat2/hg38 -1 ${0}_1_val_1.fq.gz -2 ${0}_2_val_2.fq.gz 2> ${0}.log | samtools sort -@ 2 -o ${0}.bam 

参数介绍

• -p 指定计算的线程数，可加速计算
• -x 指定索引文件（在“自建索引（hisat2）”中生成）。千万注意，最后一定要加上前述指定的索引前缀，此处为“hg38”。否则会报错。
• -1，-2 分别后接双端序列中一个 单独该命令会生成*.sam文件，但是后接了管道符“|”，故*.sam文件不会以文件形式保存下来，节省内存。
• “2>” 作用为将命令的报错重定向到后续指定的日志文件中，类似地， “1>” 作用为将正确的输出重定向到日志文件。“> a.txt 2>&1”为将所有输出都重定向到a.txt文件中。
• samtools 为将*.sam文件转化为二进制的*.bam文件，节省内存

1. 根据对软件参数的了解，发现只要通过一定的方法实现自动替换-1,-2后接的参数即可实现批量自动化的目的
2. 观察输入文件“SRR11618624_1_val_1.fq.gz SRR11618624_2_val_2.fq.gz SRR11618635_1_val_1.fq.gz SRR11618635_2_val_2.fq.gz SRR11618647_1_val_1.fq.gz SRR11618647_2_val_2.fq.gz”发现，只要提取第一个下划线的文件名前缀对命令行参数进行填充即可
3. 如下shell命令即可实现：

ls *fq.gz | xargs -t -n1 sh -c 'echo ${0%%_*} ' | sort | uniq

拿到所有文件前缀

命令介绍

• ls *fq.gz 列出所有fq.gz结尾的文件
• xargs -t -n1

• -t 在终端中将要执行的命令显示出来
• -n 每次处理几个前面传过来的文件。注：“-n1与-n 1作用相同”，后续可以用 {0}、${1}分别指定第一个、第二个文件。**

• %%_* 为字符串切分命令，作用为：对于“SRR11618624_1_val_1.fq.gz”字符串，仅保留最左边“_”左边的字符，即“SRR11618624”。其余字符串切分命令见链接。 linux中字符串截取的八种方法
• 可以想象，至此3对6个序列文件会得到6个前缀，其中两两重复。
• 经过sort(排序)及uniq(去重)之后将得到唯一的三个文件前缀“SRR11618624，SRR11618635， SRR11618647”

综上可得： 作用：批量对指定目录下所有序列文件进行比对

ls *fq.gz | xargs -t -n1 sh -c 'echo ${0%%_*} ' | sort | uniq | xargs -t -n 1 nohup sh -c 'hisat2 -p 2 -x $HOME/Project4/source_data/cleandata/refANDgtf/index_by_hisat2/hg38 -1 ${0}_1_val_1.fq.gz -2 ${0}_2_val_2.fq.gz 2> ${0}.log | samtools sort -@ 2 -o ${0}.bam ' &

注： 命令行中出现了5次${0}，第一次为指代完整的文件名，如：“SRR11618624_1_val_1.fq.gz”，其余均为指代提取到的前缀，如：“SRR11618624”

未解决的问题：

原本hisat2的“-x”参数我准备通过引用来解决，即在此命令执行之前先定义索引位置hisat_index=" {hisat_index}”或“ ”方式指定的地址识别不了。后面发现通过“ HOME/Project4/source_data/cleandata/refANDgtf/index_by_hisat2/hg38”。猜想可能是全局变量跟局部变量的关系，这个问题目前仍不理解也没解决。朋友们有知道的吗，可以在评论区教教我。

其他

• 如何得到基因组注释文件（*.gff或*.gft）与参考基因组文件（*.fa或*.fa.gz）？
• 找到一篇汇总了目前主流数据库（UCSC/ensemble/NCBI/gencode）中下载两个文件的方法的博文。 构建index所需的参考基因组以及各种版本的注释文件
• 我是通过NCBI下载两个文件，故以该数据库为例讲解：

1. 首先进入NCBI官网：https://www.ncbi.nlm.nih.gov/genome
2. 选择物种

3. 选择需要的版本，需要两个文件的版本一致

• 推文多平台同步发布，公众号内容食用更佳
• 更多内容，请关注微信公众号“生信矿工”
• 如有意见或建议可以在评论区讨论

公众号二维码

本文由 mdnice 多平台发布