搜索
查看: 4310|回复: 2

[Other] Biostar21、22-使用samtools和FreeBayes进行变异的calling并用snpEff...

[复制链接]

33

主题

46

帖子

230

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
230
发表于 2017-8-26 20:28:58 | 显示全部楼层 |阅读模式
注:本文是生信媛微信公众号原创文章
作者:István Albert
原文链接:http://mp.weixin.qq.com/s/vR6gPWiTRmkonuaN6K_Xaw

翻译:生物女博士
注:本系列课程翻译已获得原作者授权。

#为作者的注释
#*为译者的注释




Lecture 21 - 使用samtools进行变异的calling

# 建一个转换和排序的BAM文件的“快捷方式”。
alias tobam='samtools view -b - | samtools sort -o - tmp '

# 安装bcftools.
cd ~/src
tar jxvf bcftools-1.1.tar.bz2
cd bcftools-1.1
make
ln -s ~/src/bcftools-1.1/bcftools ~/bin

# 使用Lecture18的比对脚本
bash compare.sh

# 有参考基因组地Pileup
samtools mpileup -f ~/refs/852/NC.fa -Q 0 bwa.bam | more

# 生成基因型,然后把变异call出来。
# Samtools是为二倍体的人类基因组设计的,可能会有与之对应的隐含假设包含其中。
samtools mpileup -Q 0 -ugf ~/refs/852/NC.fa bwa.bam | bcftools call -vc > samtools.vcf

# 查看call出来的snp。
# snp calling是如何工作的?
# 以我写的一个DIY的call snp程序作为简单的演示(见网站)
#* 代码在本章节结尾
samtools mpileup -f ~/refs/852/NC.fa -Q 0 bwa.bam > bwa.pileup
cat bwa.pileup | python snpcaller.py > diy.vcf

# 我简单版的程序跟samtools用默认设置时得到一样的结果——感觉如何?

# 欢迎来到你的下一个文件格式。
# 按列分割文件。
cat samtools.vcf| grep -v "##" | cut -f 1,2,3,4,5 | head

# 安装Freebayes
# Mac上我们需要另一个程序来编译freebayes。
brew install cmake

# 现在来安装FreeBayes
#* 使用git clone需要先安装git。
#*也可以直接到https://github.com/ekg/freebayes.git上下载到本地,解压,然后将文件夹上传到服务器的~/src目录下进行安装。

cd ~/src
git clone --recursive git://github.com/ekg/freebayes.git
cd freebayes

# 在Mac上,下边的命令会报错。
make

# 但第二次运行就正常了,暂不清楚什么原因。程序似乎也是可以用的——但是可能有些功能是缺失的。
make
ln -sf ~/src/freebayes/bin/freebayes ~/bin

# 现在回到主目录,用FreeBayes来call SNP。
cd ~/edu/lec21/

# 程序该怎么用?
freebayes

# 用它来call snps。
freebayes -f ~/refs/852/NC.fa bwa.bam > freebayes.vcf

# 可视化结果
#*可以导入IGV查看。

DIY snp caller
#*这是一个python编写的程序,python用缩进来标明成块的代码,如果缩进有问题,程序也会出错。而微信公众号排版有时候会有一些问题,故建议查看原英文网站,确认自己的这个python脚本没问题。
[Python] 纯文本查看 复制代码
# DIY snip caller
import sys
from collections import defaultdict
# This code was written in about 20 minutes to demonstrate
# a really simple snp caller. As it is it calls only
# homozygous mutation defined as over 50% of calls.

print "##fileformat=VCFv4.2"
print "\t".join("#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT bwa.bam".split())

for line in sys.stdin:
    # Strip the newline then split by tabs.
    column = line.strip().split("\t")
    
    # The fifth base contains the basecalls# (zero based counting).
    bases = column[4]
    
    # The quals column will tell us how many bases are called.
    quals = column[5]
    
    # Upper case everything.
    # We do lose strand information with that.
    bases = bases.upper()    count = defaultdict(int)  
    for c in bases:        count[c] += 1

    # How many bases cover the site
    depth = float(len(quals))
    
    for m in "ATGC":
        if count[m]/depth > 0.50:
        # We got a homozygous mutations
        # produce the CVF records
        # Collect output into an array.

        info = "DP=%s" % (int(depth))        format = "GT:PL"
        values = "1/1:30"    
        out = [column[0], column[1], ".", column[2], m, 30, ".", info, format, values]
        
        # Make all elements of the array a string.
        out = map(str, out)
        
        # Join the fields of the output array by tabs
        # Print the joined fields.
        print "\t".join(out)


Lecture 22 - 预测变异的影响
# 安装、创建以及链接bioawk
cd ~/src

#* 使用git clone需要先安装git。
#* 也可以直接到https://github.com/lh3/bioawk.git上下载到本地,解压,然后将文件夹上传到服务器的~/src目录下进行安装。

cd bioawk
#* 译者下载后的名字为bioawk-master,在此步骤将命令对应更改为cd bioawk-master/ 即可,其他一样。
make
ln -sf ~/src/bioawk/bioawk ~/bin

# bioawk有自带的一个手册。
man ~/src/bioawk/awk.1

# Bioawk知道怎样处理生物信息类型数据
# 它不仅能按列分割数据,还可以识别列名
# 通过运行compare.sh脚本生成模拟reads并比对后(详情参见之前的课程),现在用bioawk来处理它们。
bioawk -c help

cat r1.fq | bioawk -c fastx ' { print $seq } ' | head -1
cat mutations.gff | bioawk -c gff ' { print $seqname, $start } ' | head -1

# 跟下边写法比起来,bioawk使得代码变得更有可读性:
cat mutations.gff | grep -v "#" | awk -F '\t' -v OFS='\t' ' { print $1, $5-$4+1 } ' | head -1

# 我们可以这么写:
cat mutations.gff | bioawk -c gff ' { print $seqname, $end-$start+1 } ' | head -1

# 从现在开始,我们将使用bioawk。我们迟点将会涉及seqtk、tabtk、tabix。
# 我们来下载snpEff
cd ~/src
unzip snpEff_latest_core.zip

# 跟trimmomatic和readseq一样,我们可以写个脚本来运行snpEff。
# 简便起见,我们先建一个别名(alias)。
alias snpEff='java -jar ~/src/snpEff/snpEff.jar'

# 运行snpEff
snpEff

# 你可能会看到一个报错Y:UnsupportedClassVersionError
# 这意味着你需要把你的java程序升级到至少是版本7。
java -version

# 搜索下载以及安装Java JDK (Java Development Kit, 而不是JRE, Java Runtime Environment)

# 现在snpEff应该可以正常运行了。
snpEff

# 有一些预建好的snpEff数据库
snpEff databases | more

# 搜索ebola的数据
snpEff databases | grep ebola

# 下载ebola的数据库
# -v 命令使得过程变得详细:打印了更多的信息。
snpEff download -v ebola_zaire

# 但是,这些注释是什么?
snpEff dump ebola_zaire

# 好,所以注释是跟另一个基因组相关的,而不是我们现在用的这个。我们需要获取那个基因组。

# 同样获取genbank文件,并提取编码序列到gff文件。
readseq -format=FASTA -o ~/refs/852/KJ660346.fa KJ660346.gb
readseq -format=GFF  -o KJ660346.gff KJ660346.gb

# query为NC.fa里的1999年的旧基因组。
# 通过bwa和samtools对基因组进行建立索引。
bash align-genomes.sh ~/refs/852/KJ660346.fa ~/refs/852/NC.fa

# 现在我们需要创建一个脚本,把两个文件作为输入并生成一个VCF文件。
# 运行snpEff
java -jar ~/src/snpEff/snpEff.jar -v  ebola_zaire aln.vcf > annotated.vcf


Align two genomes
[Shell] 纯文本查看 复制代码
# 通用的单端比对软件和snp caller。
# 第一个参数是参考序列,第二个参数是query

# 这可以在没有参数传递时终止程序。
set -ue

# 从命令行传来参数
# 第一个参数是第一个参考基因组文件
GENOME1=$1

# 第二个参数是第二个参考基因组
GENOME2=$2
#*bash align-genomes.sh ~/refs/852/KJ660346.fa ~/refs/852/NC.fa
#*上边这行代码中 ~/refs/852/KJ660346.fa便是传入到脚本align-genomes.sh的第一个参数,即$1,而~/refs/852/NC.fa则是第二个传入的参数$2

# 这只需运行一次
#*如果之前已经index(索引)好了,可以不做这一步
bwa index $GENOME1

# 运行bwa并生成bwa.sam文件
bwa mem $GENOME1 $GENOME2 > aln.sam

# 假如你已经安装了lastz比对软件,你也可以用它。
#lastz $GENOME1[unmask] $READS[unmask] --format=sam > aln.sam

samtools view -Sb aln.sam > tmp.bam
samtools sort -f tmp.bam aln.bamsamtools index aln.bam

# 去除中间文件。
rm -f tmp.bam

# Call snps
samtools mpileup -f $GENOME1 -ug aln.bam | bcftools call -vm > results.vcf

欢迎到微信公众号订阅我们
生信媛
bio_sxy





上一篇:Biostar:课程19、20-利用samtools mpileup和bcftools进行SNP calling
下一篇:bam文件里面的测序质量值问题 phred64 vs phred33
回复

使用道具 举报

0

主题

7

帖子

51

积分

注册会员

Rank: 2

积分
51
发表于 2017-10-13 18:02:31 | 显示全部楼层
做出来 diy.vcf;samtools.vcf;freebayes.vcf的结果好像不一样,不知道什么原因?特别是diy.vcf和samtools.vcf之间按照原文理解应该是一样的,也有些许差异,freebayes的差异就更大了。。。。

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

33

主题

46

帖子

230

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
230
 楼主| 发表于 2017-10-21 10:33:09 | 显示全部楼层
二口 发表于 2017-10-13 18:02
做出来 diy.vcf;samtools.vcf;freebayes.vcf的结果好像不一样,不知道什么原因?特别是diy.vcf和samtools.v ...

不同软件,因为算法的不同,call出来的SNP是不一样的。
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-7-24 09:15 , Processed in 0.037172 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.