linux序列比对程序,序列比对工具的对比

# 下载并安装比对软件bowtie2cd ~/src

# Mac OSX上用:curl -OL http://downloads.sourceforge.net/project/bowtie-bio/bowtie2/2.2.4/bowtie2-2.2.4-macos-x86_64.zip

# Linux上用:#curl -OL http://downloads.sourceforge.net/project/bowtie-bio/bowtie2/2.2.4/bowtie2-2.2.4-linux-x86_64.zip

# 这个已经是个二进制文件了,所以可以直接执行。unzip bowtie2-2.2.4-macos-x86_64.zip

# 把比对软件以及相关程序链接到bin文件夹。ln -s ~/src/bowtie2-2.2.4/bowtie2 ~/bin/

ln -s ~/src/bowtie2-2.2.4/bowtie2-align ~/bin/

ln -s ~/src/bowtie2-2.2.4/bowtie2-build ~/bin/

# 对参考基因组建立索引。

# 这是bowtie2的索引。bowtie2-build ~/refs/852/NC.fa NC.fa

# 把突变弄成一定的格式,使得可以在浏览器上展示。

# 把它变为gff文件。

#*mutations.txt这个文件是在lecture15中生成的。并且存放在目录~/edu/lec15下。cat mutations.txt | awk ' {print $1, "wgsim", "mutation", $2, $2, ".", "+",  ".", "." }' > mutations.gff

# 运行比对程序。

#*脚本见文末bash compare.sh

# 修改脚本使之生成有较大的错误率的fastq文件。

# 当你修改了错误率后,计算一下比对上的reads数目。samtools view -cF 4 bwa.bam

samtools view -cF 4 bow.bam

最终脚本compare.sh# 对比两个比对软件的输出。

# 使用方法: bash compare.sh

# 在出错的地方停止运行。

set -ue

# 数据文件名。

READ1=r1.fq

READ2=r2.fq

# 参考序列文件。两个比对软件需要分别对它进行索引。

REFS=~/refs/852/NC.fa

# 从基因组中生成模拟reads。

# 编辑错误率并重新运行这个pipeline。

wgsim -N 10000 -e 0.1 $REFS $READ1 $READ2 > mutations.txt

# 将mutations文件转成GFF文件。

cat mutations.txt | awk -F '\t' ' BEGIN { OFS="\t"; print "##gff-version 2" } { print $1, "wgsim", "mutation", $2, $2, ".", "+",  ".", "name " $3 "/" $4 }' > mutations.gff

# 我们需要添加一个read组到mapping中。

GROUP='@RG\tID:123\tSM:liver\tLB:library1'

# 运行bwa并创建bwa.sam文件。

bwa mem -R $GROUP $REFS $READ1 $READ2 > bwa.sam

# 运行bowtie2并创建bow.sam文件。

bowtie2 -x $REFS -1 $READ1 -2 $READ2 > bow.sam

# 调整bowtie2# bowtie2 -D 20 -R 3 -N 1 -L 20 -i S,1,0.50 -x $REFS -1 $READ1 -2 $READ2 > bow.sam

# 对每个sam文件,都转换成bam(sam的二进制)格式。

for name in *.sam;

do

samtools view -Sb $name > tmp.bam

samtools sort -f tmp.bam $name.bam

done

# 删除中间文件。

rm -f tmp.sam tmp.bam

# 修改名字,使得他们不再含有两个后缀。

mv bwa.sam.bam bwa.bam

mv bow.sam.bam bow.bam

# 对数据进行索引。

for name in *.bam

do

samtools index $name

done

# 删除这个程序生成的所有文件。

# rm -f *.bam *.bai *fq *.txt *.sam *.gff

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值