J

Jay2Echo

V1

2022/12/11阅读:18主题:默认主题

[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. 选择需要的版本,需要两个文件的版本一致 NCBI网站下载方法介绍

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

分类:

前端

标签:

工具介绍

作者介绍

J
Jay2Echo
V1