背景
前面介绍了如何从sra文件得到fastq文件,并利用fastqc、multiqc生成质控报告以及使用trim_galore去除低质量碱基
[Linux|生信]project4_01:批量下载sra文件并转化为fastq文件
[Linux|生信]project4_02:质控过滤
今天介绍一下如何自建序列索引以及序列比对
Just Do It!
输入文件准备
-
基因组注释文件(*.gff或*.gft) -
参考基因组文件(*.fa或*.fa.gz) -
前述经过质控的序列文件( .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
![生成的索引文件](https://i-blog.csdnimg.cn/blog_migrate/cdaa4abe5b85c4aeb774ec51c0c7f28b.png)
参数说明
extract_exons.py与extract_splice_sites.py此两个Python脚本为hisat2软件内置的脚本,成功安装即可使用 参考基因组与基因组注释文件均位于当前目录的上一级目录,操作时注意修改相应文件路径 hisat2-build用法介绍
-p: 指定计算的线程数,可加速运算 --ss 和 --exon 后接文件为前述生成 .fna为参考基因组文件,正常为.fa。我下载的为*.fna格式,打开看了一下,里面存在大量的不确定碱基被表注为“N”,猜想*.fna文件为包含不确定碱基“N”的参考基因组文件。不影响后续操作。 最后的hg38为自定义的索引前缀,运行完发现生成了8个索引文件
比对
![质控后得到的三对双端序列文件](https://i-blog.csdnimg.cn/blog_migrate/a74d15c6a47814179bd3b699d7384c39.png)
-
首先将质控后的双端序列文件拷贝至一个空文件夹,以确保其他*.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,-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”发现,只要提取第一个下划线的文件名前缀对命令行参数进行填充即可 -
如下shell命令即可实现:
ls *fq.gz | xargs -t -n1 sh -c 'echo ${0%%_*} ' | sort | uniq
![拿到所有文件前缀](https://i-blog.csdnimg.cn/blog_migrate/6236cfa89ce088a5e781b629edc994c7.png)
命令介绍
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下载两个文件,故以该数据库为例讲解:
-
首先进入NCBI官网:https://www.ncbi.nlm.nih.gov/genome -
选择物种
3. 选择需要的版本,需要两个文件的版本一致
-
推文多平台同步发布,公众号内容食用更佳 -
更多内容,请关注微信公众号“生信矿工” -
如有意见或建议可以在评论区讨论
![公众号二维码](https://i-blog.csdnimg.cn/blog_migrate/eca9947413334222324ea12a48ca86f8.png)
本文由 mdnice 多平台发布