# 下载并安装比对软件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