shell脚本实现基因重测序和变异位点检测(附源代码)(来源于知乎:我爱小徐子)

0 人赞了该文章

基因组重测序和变异的这一节,之前有花时间学习过,所以整个过程,心中也大致有一些脉络,下面就当是做一个回顾。 下面引用一下,谈老师ppt里面的一张图:

这张图就是整个重测序的流程。


之前在学习完毕之后,有自己写过一个shell脚本,能够将整个重测序的流程走一遍,并输出正确的结果,当时老师提意见让我把它写成一个更加专业的脚本,把help文档之类的都写上去,,本人比较懒惰,最近又到了复习备考的季节,所以就没有继续做下去,不过可以看看实现的效果,如下图: 首先按照usage输入需要的文件,脚本自动检测是否安装了相关软件

继续运行

进入一个小的循环系统,可以选择一些相应的功能

最终得到最后的结果



上面的图片就是整体shell脚本运行的效果,是交互式的,可以选择一些不同的功能。

源代码在文末


下面我来大致理一下整个重测序变异位点检测的思路: 从文件格式的角度来考虑,整个过程无非就是从: .fq(需要以.fna格式文件做参考基因组)->.sam->.bam->.sorted(.bam的排序文件)->.bcf(二进制文件)->vcf(文本格式)->.gff(自定义数据库时需要使用的参考文件)->.variant_function/.exonic_variant_function文件


从需要的软件来考虑: bwa->samtools->bcftols->GenePred->annovar


大致的思路过程就是上述,整体上不难,但是有一个小难点就是自定义数据库的地方,下面详细讲解一下: vcf文件是用于描述SNP,InDel和SV结果的文本文件。VCF格式在GATK软件中得到很好的支持。在生成.vcf文件之后,我们可以使用bcftools对它进行过滤,得到.snps.vcf文件。 下面我们需要做的是对变异基因进行注释,使用的软件是Annovar。 annovar的下载地址:

http://www.openbioinformatics.org/annovar/download/0wg xR2rIVP/annovar.latest.tar.gz

这里再次引用谈老师的一张ppt:


我们在这里需要用到的pl脚本是:

  1. convert2annovar.pl,它的作用主要是将.vcf文件格式转化为.avinput格式,用于程序的后续处理;
  2. annotate_variation.pl,它用来注释,使用她的时候注意格式:usage: annotate_variation.pl --geneanno --dbtype refGene -buildver 7942-genome-new test_bwa_7942.snps.avinput ~/Biosofts/annovar/humandb/
    特别注意,~/Biosofts/annovar/humandb/是自定义数据库所在的位置 生成最终的结果,结果的格式为:.avinput.variant_function文件和 .avinput.exonic_variant_function文件

3.retrieve_seq_from_fasta.pl 用于自行建立其他物种的转录本

之后就是使用gff文件自定义数据库:

GFF全称为general feature format,是Sanger研究所 定义的,一种简单的、方便的对于DNA、RNA以 及蛋白质序列的特征进行描述的数据格式。主要是用来注释基因组,告诉我们序列的那个区 域是基因、RNA等。

./gff3ToGenePred将.gff格式的转化为.txt文件,之后使用retrieve_seq_from_fasta.pl脚本生成参考数据库,最终所有的中间结果加在一起,得到的数据如图:



以上就是整个重测序变异检测的大致过程。 下面给出我的shell源代码(代码里面有注释):

#!/bin/bash

:<<!
常用重测序数据分析流程shell脚本实现整个过程

##大致思路:
1.需要的软件
bwa
samtools
bcftools(属于samtools)
GenePred
annovar
2.步骤
1)下载参考序列.fna,和需要比对的fastq格式的reads,(.fq)
2)获得.sam文件(使用bwa)
3)samtools转化.sam为.bam
4)对.bam进行sorted整理
5)进行功能选择:1,2,3,4
6)生成.bcf文件
7)使用bcftools检测基因变异位点
8)变异位点过滤
9)使用annovar
!

isInstalled ()
{
  indelString=""
  echo -e  "正在检测您是否安装了必须软件..."
  if [  -z  `which bwa` ]
  then
  echo "您未安装bwa软件,或者未将其加入环境变量"
  else
   echo "bwa检测成功..."
  fi
  if [ -z  `which samtools` ]
  then
  echo "您未安装samtools软件,或者未将其加入环境变量"
  else
  echo "samtools检测成功..."
  fi
  if [  -z  `which bcftools` ]
  then
  echo "您未安装bcftools软件,或者未将其加入环境变量"
  else
  echo "bcftools检测成功..."
  fi
  if [  -z   `which gff3ToGenePred` ]
  then
  echo "您未安装gff3ToGenePred软件,或者未将其加入环境变量"
  else
  echo "gff3ToGenePred检测成功..."
  fi
  if [  -z  `which annotate_variation.pl` ]
  then
  echo "您未安装annovar软件,或者未将annovar文件夹加入环境变量"
  else
  echo  "annovar检测成功..."
  fi
}
#isInstalled
#得到参考基因序列,需要3个参数
fna=$1
fq1=$2
fq2=$3
gff=$4

getSamSeq()
{
  bwa  index $fna
  bwa mem $fna $fq1 $fq2 > ./SamSeq.sam
  #echo "$1,$2,$3,$4"
}

#getSamSeq
#得到bam序列,需要  1个参数
getBamSeq()
{
  samtools  faidx $fna
  a=` find ./ -name "*.fai"`
  echo "$a"
  samtools  view -bhS -t $a  -o ./BamSeq.bam  ./SamSeq.sam
  samtools  sort  ./BamSeq.bam -o  ./BamSeq.bam.sorted
  samtools  index  ./BamSeq.bam.sorted
}

#getBamSeq
#选择功能

selectFun()
{
  echo -e "请输入您选择的功能的序号:\n1.显示reads比对基因组的情况\n2.查看reads  mapping到基因组的情况\n3.测试>参考基因组每个位点或一段区域的测序深度\n4.统计比对的结果\n5.退出系统\n"
  read -p "请输入序号:"  num
  case "$num"  in
  "1")
 samtools  tview  ./BamSeq.bam.sorted  $fna
  ;;
  "2")
 echo -e "自行下载Artemis在windows系统下查看\n"
  ;;
  "3")
  samtools  depth  ./BamSeq.bam.sorted >depth.txt
  echo  -e "以将结果输出至depth.txt文件,请退出查看\n"
  ;;
  "4")
   samtools  flagstat  ./BamSeq.bam.sorted
  ;;
  "5")
  exit
  ;;
  *)
echo "您的输入有误"
  ;;
esac

}

#生成bcf文件
getBcfFile()
{
  bcftools mpileup -f $fna -O "u"  ./BamSeq.bam.sorted   -o  ./Seq.bcf
  bcftools call  -vm ./Seq.bcf -o  ./Seq.variants.bcf
  echo "进行snp与indel变异位点检测..."
  bcftools view -v  snps,indels  ./Seq.variants.bcf >./Seq.snps.vcf
  echo "对变异位点进行过滤:变异质量大于20,深度大于5"
 bcftools filter  -o ./Seq.snps.filtered.vcf  -i 'QUAL>20&&DP>5' ./Seq.snps.vcf

}
#getBcfFile
#变异基因检测
getAnnovar()
{
  convert2annovar.pl -format  vcf4  ./Seq.snps.filtered.vcf >./Seq.avinput
  annotate_variation.pl --geneanno --dbtype  refGene  --buildver  Seq  Seq.avinput  ./

}
#getAnnovar
#自定义数据库
myBioDb()
{
  gff3ToGenePred -useName  $gff  Seq_refGene.txt
  retrieve_seq_from_fasta.pl -format refGene  -seqfile  $fna -outfile Seq_refGeneMrna.fa  ./Seq_refGene.txt
}
#myBioDb
:<<!
函数定义结束,下面开始主代码,
整个脚本需要输入4个参数,$1是.fna文件  $2,$3是fq.gz文件  $4是.gff文件

!
main()
{
  isInstalled
  getSamSeq
  getBamSeq
  getBcfFile
  myBioDb
  getAnnovar
    echo "您可以选择一些功能:"
  while :
  do
  selectFun

  done
}
main


更多优质内容,请关注微信公众号:生物信息与python

欢迎各位生信大牛来访生物信息社区:biocoder.cn(功能还在完善中,最近要考试了,寒假更新)

django学习小组:921464927


评论

暂无评论