序列比对工具的对比

2017/01/25来源：生信媛微信公众号评论6,444

# 下载并安装比对软件bowtie2

cd ~/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

发表评论