Biostar：课程15、16
翻译：生物女博士
注：本系列课程翻译已获得原作者授权。
#为作者的注释
#*为译者的注释
Lecture 15 - BAM文件和samtools
# 安装一个短read模拟器。Heng Li写的wgsim。
cd ~/src
git clone https://github.com/lh3/wgsim.git
cd wgsim
gcc -g -O2 -Wall -o wgsim wgsim.c -lz -lm
ln -s ~/src/wgsim/wgsim ~/bin/wgsim
# 检查一下是否可以用。
wgsim
# 现在，下载和安装samtools。
# 把samtools链接到~/bin目录。
cd ~/src
curl -OL http://sourceforge.net/projects/samtools/files/samtools/1.1/samtools-1.1.tar.bz2
tar jxvf samtools-1.1.tar.bz2
cd samtools-1.1
make
ln -s ~/src/samtools-1.1/samtools ~/bin/
# 检查一下是否可以用。
samtools
# 查看samtools的手册。
man ~/src/samtools-1.1/samtools.1
# 我们将会利用一个突变的参考基因组来生成一些短reads，并map（回贴？）回去。
mkdir ~/edu/lec15
cd ~/edu/lec15
# 生成一些数据，并把输出保存到一个文件。
wgsim -N 10000 ~/refs/852/ebola-1999.fa read1.fq read2.fq > mutations.txt
# 运行比对。
bwa mem ~/refs/852/ebola-1999.fa read1.fq read2.fq > results.sam
# 开始转换成BAM文件。
#*运行的时候请注意，这里用的samtools版本比较旧，新版本有的命令已经改变了。
samtools view -Sb results.sam > temp.bam
# 对比对结果进行排序。
samtools sort -f temp.bam results.bam
# 索引比对结果。
samtools index results.bam
#*做完了这一步，如果想方便查看一下比对结果，可以下载一个IGV软件，并将参考序列、bam文件和bam文件的索引后生成的一个文件导入IGV即可。
# 通常，我们会把所有的命令放到一个脚本中并运行。
# 查看文末的脚本。
#*原文并未将这个脚本放上来。但是后边的课程中会有更加完善的脚本。
#*本行命令暂时不运行：bash align.sh read1.fq read2.fq results.bam
# 使用samtools过滤。
# -f 匹配标识（flag）（保留匹配的标识对应的reads）
# -F 过滤标识（flag）（去除匹配标识对应的reads，保留剩余的）
#*在上两个lecture中有说过sam格式，sam每一列有不同含义，其中一列（第二列）叫标识（flag），不同数字有不同含义
#*下边内容引用自参考资料1
FLAG, 概括出一个合适的标记，各个数字分别代表
1 序列是一对序列中的一个
2 比对结果是一个pair-end比对的末端
4 没有找到位点
8 这个序列是pair中的一个但是没有找到位点
16 在这个比对上的位点，序列与参考序列反向互补
32 这个序列在pair-end中的的mate序列与参考序列反响互补
64 序列是 mate 1
128 序列是 mate 2
假如说标记为以上列举出的数目，就可以直接推断出匹配的情况。假如说标记不是以上列举出的数字，比如说83=（64+16+2+1），就是这几种情况值和。
# 比对到反向互补链的reads。
samtools view -f 16 results.bam
# 比对到正向链的reads。
samtools view -F 16 results.bam
# -c 可以用来计数
samtools view -F 16 results.bam | wc -l
samtools view -c -F 16 results.bam
# 用质量来过滤。BWA 的mapping质量为0时，表示reads是map到多个位置。而 q>=1表示该reads是单一mapping。
samtools view -c -q 1 results.bam
# 统计高质量的比对数目。
samtools view -c -q 40 results.bam
# 参考基因组的名字为 gi|10313991|ref|NC_002549.1|
# 总是手动输入会有些乏味，我们可以创建一个简称。
CHR='gi|10313991|ref|NC_002549.1|'
# 截取部分数据。
samtools view results.bam $CHR:1-100
# 查看文件的特定某列
samtools view results.bam $CHR:1-100 | cut -f 4 | more
Lecture 16 - Converting data with ReadSeq
# 安装文件格式转换器。
# 为之创建一个目录。
mkdir -p x
cd ~/src/readseq
curl -OL http://iubio.bio.indiana.edu/soft/molbio/readseq/java/readseq.jar
# 跟调用trimmomatic很相似
# 将调用trimmomatic脚本修改一下。
#*在lecture8时，我们创建了一个文件叫trimmomatic来简化trimmomatic的命令执行。详情参见Biostar：课程7、8
cp ~/bin/trimmomatic ~/bin/readseq
# 用下边程序来替换：
# java -jar ~/src/readseq/readseq.jar $@
# 你也可以加上：
# java -cp ~/src/readseq/readseq.jar app $@
# 用图形用户界面来运行。
# 获取完整的GenBank记录的1999年的埃博拉基因组。
# http://www.ncbi.nlm.nih.gov/genome/4887
# 获取完整genbank文件。
mkdir -p ~/edu/lec16
cd -p ~/edu/lec16
# 获取完整数据。
#*efetch貌似有问题，可以用网页直接下载。详情参见Biostar：课程3、4
efetch -db nucleotide -id NC_002549.1 -format gb > NC.gb
# 提取数据未GFF（Generic Feature Format）格式。
# 自动获取输入文件的格式。
readseq -format=GFF NC.gb
# 你也可以设置输出文件的名字。
readseq -format=GFF -o NC-all.gff NC.gb
# 提取为fasta文件。
readseq -format=FASTA -o NC.fa NC.gb
# 从GFF文件中提取“gene”的行
#* \t表示的是tab分隔符
cat NC-all.gff | egrep '\tgene\t'
# 同时，前两行也是我们要的。
echo '##gff-version 2' > NC-genes.gff
cat NC-all.gff | egrep '\tgene\t' >> NC-genes.gff
# 染色体坐标不匹配！我们需要去把它修复一下。
# 我们需要重做比对，或者交换现有比对或特征的坐标。
# 重新索引和比对数据。
cp NC.fa ~/refs/852/
bwa index ~/refs/852/NC.fa
cp ../lec15/*.fq .
# 编辑align.sh使其能用于新文件，并重新运行。
bash align.sh read1.fq read2.fq results.bam
# 把数据导入IGV并查看。
# 查看100-500bp区域的深度。
# 输出为序列号、基因组索引和覆盖度。
samtools depth -r NC_002549:100-500 results.bam
参考资料1：http://www.cnblogs.com
Biostar非常适合有生物学基础且想自学生信的这部分人入门。
Biostar课程文章系列：
【课程预告】手把手教你入门生信——The Biostar Handbook
Biostar：课程1、2
欢迎关注我们