本地blast使用
- 下载及安装：
curl -O ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.6.0+-x64-linux.tar.gz
tar zxvf ncbi-blast-2.6.0+-x64-linux.tar.gz
#添加路径
echo 'export PATH=$PATH:~/src/ncbi-blast-2.6.0+/bin' >> ~/.bashrc
source ~/.bashrc
- 使用：
blast与大多数序列比对软件一样，需要对参考序列创建一个哈希结构的数据库：
#参数：
# -dbtype <String, `nucl', `prot'>
# Molecule type of target db
# -in <File_In>
# Input file/database name
# Default = `-'
makeblastdb -in ~/refs/KM233118.fa -dbtype nucl
序列比对
# 我们可以修改输出格式。 当碱基对碱基的比对信息不需要时，我们可以采用表格格式。
blastn -db refs/KM233118.fa -query query.fa -outfmt 6
# 你也可以自己定义输出格式。
blastn -db refs/KM233118.fa -query query.fa -outfmt '6 qseqid sseqid qlen length mismatch'
# blast有不同的搜索策略或task。
#*task里的参数有'blastn'，'blastn-short' ，'dc-megablast'，'megablast'， 'rmblastn'，默认的是'megablast'
#*使用不同的task，出来的结果内会有文献引用，有兴趣的朋友可以去看一下。
#*blastn：需要11bp以上的精确匹配。
#*blast-short：对于小于50bp的序列最优。
#*megablast：一般用于非常相似的序列。
#*dc-megablast:dc（discontiguous，不连续的），一般用于找距离远（相似性低的，比如种间同源）的序列。
blast可以对多条序列进行比对
# 你在目标数据库中，可以有多条序列。
esearch -db nucleotide -query PRJNA257197 | efetch -db nucleotide -format fasta > refs/all-genomes.fa
# 对所有的基因组创建blast数据库。
makeblastdb -in refs/all-genomes.fa -dbtype nucl
# 随意挑取基因组的任意区域。
efetch -db nucleotide -id KM233118 -format fasta -seq_start 1 -seq_stop 1000 > 1000.fa
#*大家可以试一下
#*blastn -db refs/all-genomes.fa -query 1000.fa -outfmt 6
blastdbcmd, blastp, blastx, tblastn 的用法
# blastdbcmd 有很多参数，可以对Blast数据库内容进行查询和格式设定。
# 列出数据库的所有内容。
blastdbcmd -db refs/NC_002549.fa -entry 'all' | head -3
# 获取数据库的特定条目。
blastdbcmd -db refs/NC_002549.fa -entry '10313991' | head -3
#*报错，暂未解决。
# 获取数据库中的一定范围的核苷酸条目。
blastdbcmd -db refs/NC_002549.fa -entry 'all' -outfmt '%a %s' -range 1-10 -strand minus
#*%a表示accession（访问号），%s 表示 sequence data（序列数据）
# 设定数据库内容的格式。
blastdbcmd -db refs/NC_002549.fa -entry 'all' -outfmt '%a %l'
#*%l 表示 sequence length（序列长度）
# 列出每个蛋白及其长度。
blastdbcmd -db refs/NC_002549-prot.fa -entry 'all' -outfmt '%a %l'
# 从蛋白数据库中提取一个条目并存到一个文件中。
blastdbcmd -db refs/NC_002549-prot.fa -entry 'NP_066243.1' > query-p.fa
blastp是蛋白序列比对蛋白数据库，与blastn用法类似，也是需要先进行数据库的构建，再进行比对：
blastp -db refs/NC_002549-prot.fa -query query-p.fa
tblastn是将蛋白序列比对到核苷酸数据库。
efetch -db protein -id NP_066243.1 -format fasta > NP_066243.fa
tblastn -db refs/NC_002549.fa -query NP_066243.fa | more
blastx是核苷酸序列比对到蛋白数据库
blastx -db refs/NC_002549-prot.fa -query nucleotide.fa | more
BWA(Burrows-Wheeler Aligner)
BWA主要是将reads比对到大型基因组上，主要功能是：序列比对。首先通过BWT(Burrows-Wheeler Transformation，BWT压缩算法)为大型参考基因组建立索引，然后将reads比对到基因组。特点是快速、准确、省内存。由三种类似算法组成：BWA-backtrack，BWA-SW和BWA-MEM。首推BWA-MEM。
三种算法的适用范围
- BWA-backtrack：reads长度<70bp时，推荐本算法，建议输入reads长度 < 100bp。
- BWA-SW：在reads具有频繁的gap时，比对更敏感，推荐本算法。reads长度一般为70bp-1Mbp，支持long-reads，split alignment。
- BWA-MEM(首推)：在reads长度在70bp-1Mbp范围时，推荐本算法(除了上面两种情况)。支持long-reads，split alignment。
- 安装及下载
cd ~/src
curl -OL http://downloads.sourceforge.net/project/bio-bwa/bwa-0.7.15.tar.bz2
tar jxvf bwa.0.7.15.tar.bz2 #bzip2压缩，需要加-j 解压
cd bwa-0.7.15/
make
#建立软链接
ln -s ~/src/bwa-0.7.15/bwa ~/bin
- 使用
# 对基因组创建索引，这样bwa就可以进行基因组搜索。
bwa index ~/refs/852/ebola-1999.fa
# 获取埃博拉基因组序列的第一行。
head -2 ~/refs/852/ebola-1999.fa > query.fa
# 用bwa-mem命令将其map回它的基因组。
bwa mem ~/refs/852/ebola-1999.fa query.fa > results.sam
SAM文件格式
主要细节参照该网址：https://www.cnblogs.com/zdwu/p/6801063.html
# 有两种输出流。它们被称为：标准输出（standard out）和标准错误输出（standard error）。
# 下边的代码是告诉你，你可以如何重定向这两个输出到文件中。
bwa mem ~/refs/852/ebola-1999.fa query.fa 1> results.sam 2> errors.txt
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
- 下载和安装samtools。
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/
- 使用
# 生成一些数据，并把输出保存到一个文件。
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 view -Sb results.sam > temp.bam
# 对比对结果进行排序。
samtools sort -f temp.bam results.bam
# 使用samtools过滤。
# -f 匹配标识（flag）（保留匹配的标识对应的reads）
# -F 过滤标识（flag）（去除匹配标识对应的reads，保留剩余的）
#*在上两个lecture中有说过sam格式，sam每一列有不同含义，其中一列（第二列）叫标识（flag），不同数字有不同含义
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
ReadSeq:格式转换软件
- 下载与安装
cd ~/src/readseq
curl -OL http://iubio.bio.indiana.edu/soft/molbio/readseq/java/readseq.jar
- 使用
# 提取数据未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
