本地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與大多數(shù)序列比對軟件一樣疹味，需要對參考序列創(chuàng)建一個哈希結構的數(shù)據(jù)庫：
#參數(shù)：
# -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里的參數(shù)有'blastn'叙谨，'blastn-short' 许溅，'dc-megablast'瓤鼻，'megablast'， 'rmblastn'贤重，默認的是'megablast'
#*使用不同的task茬祷，出來的結果內(nèi)會有文獻引用，有興趣的朋友可以去看一下并蝗。
#*blastn：需要11bp以上的精確匹配祭犯。
#*blast-short：對于小于50bp的序列最優(yōu)。
#*megablast：一般用于非常相似的序列滚停。
#*dc-megablast:dc（discontiguous沃粗，不連續(xù)的），一般用于找距離遠（相似性低的铐刘，比如種間同源）的序列陪每。
blast可以對多條序列進行比對
# 你在目標數(shù)據(jù)庫中影晓，可以有多條序列镰吵。
esearch -db nucleotide -query PRJNA257197 | efetch -db nucleotide -format fasta > refs/all-genomes.fa
# 對所有的基因組創(chuàng)建blast數(shù)據(jù)庫。
makeblastdb -in refs/all-genomes.fa -dbtype nucl
# 隨意挑取基因組的任意區(qū)域挂签。
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 有很多參數(shù)疤祭，可以對Blast數(shù)據(jù)庫內(nèi)容進行查詢和格式設定。
# 列出數(shù)據(jù)庫的所有內(nèi)容饵婆。
blastdbcmd -db refs/NC_002549.fa -entry 'all' | head -3
# 獲取數(shù)據(jù)庫的特定條目勺馆。
blastdbcmd -db refs/NC_002549.fa -entry '10313991' | head -3
#*報錯，暫未解決侨核。
# 獲取數(shù)據(jù)庫中的一定范圍的核苷酸條目草穆。
blastdbcmd -db refs/NC_002549.fa -entry 'all' -outfmt '%a %s' -range 1-10 -strand minus
#*%a表示accession（訪問號），%s 表示 sequence data（序列數(shù)據(jù)）
# 設定數(shù)據(jù)庫內(nèi)容的格式搓译。
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'
# 從蛋白數(shù)據(jù)庫中提取一個條目并存到一個文件中。
blastdbcmd -db refs/NC_002549-prot.fa -entry 'NP_066243.1' > query-p.fa
blastp是蛋白序列比對蛋白數(shù)據(jù)庫些己，與blastn用法類似豌鸡，也是需要先進行數(shù)據(jù)庫的構建，再進行比對：
blastp -db refs/NC_002549-prot.fa -query query-p.fa
tblastn是將蛋白序列比對到核苷酸數(shù)據(jù)庫段标。
efetch -db protein -id NP_066243.1 -format fasta > NP_066243.fa
tblastn -db refs/NC_002549.fa -query NP_066243.fa | more
blastx是核苷酸序列比對到蛋白數(shù)據(jù)庫
blastx -db refs/NC_002549-prot.fa -query nucleotide.fa | more
BWA(Burrows-Wheeler Aligner)
BWA主要是將reads比對到大型基因組上涯冠，主要功能是：序列比對。首先通過BWT(Burrows-Wheeler Transformation逼庞，BWT壓縮算法)為大型參考基因組建立索引蛇更，然后將reads比對到基因組。特點是快速、準確派任、省內(nèi)存共耍。由三種類似算法組成：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
- 使用
# 對基因組創(chuàng)建索引弃甥，這樣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文件格式
主要細節(jié)參照該網(wǎng)址：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/
- 使用
# 生成一些數(shù)據(jù)搜囱，并把輸出保存到一個文件丑瞧。
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
# 開始轉(zhuǎn)換成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），不同數(shù)字有不同含義
FLAG, 概括出一個合適的標記获雕，各個數(shù)字分別代表
1 序列是一對序列中的一個
2 比對結果是一個pair-end比對的末端
4 沒有找到位點
8 這個序列是pair中的一個但是沒有找到位點
16 在這個比對上的位點薄腻，序列與參考序列反向互補
32 這個序列在pair-end中的的mate序列與參考序列反響互補
64 序列是 mate 1
128 序列是 mate 2
假如說標記為以上列舉出的數(shù)目，就可以直接推斷出匹配的情況届案。假如說標記不是以上列舉出的數(shù)字庵楷，比如說83=（64+16+2+1），就是這幾種情況值和楣颠。
# 比對到反向互補鏈的reads尽纽。
samtools view -f 16 results.bam
# 比對到正向鏈的reads。
samtools view -F 16 results.bam
# -c 可以用來計數(shù)
samtools view -F 16 results.bam | wc -l
samtools view -c -F 16 results.bam
# 用質(zhì)量來過濾童漩。BWA 的mapping質(zhì)量為0時弄贿，表示reads是map到多個位置。而 q>=1表示該reads是單一mapping矫膨。
samtools view -c -q 1 results.bam
# 統(tǒng)計高質(zhì)量的比對數(shù)目差凹。
samtools view -c -q 40 results.bam
ReadSeq:格式轉(zhuǎn)換軟件
- 下載與安裝
cd ~/src/readseq
curl -OL http://iubio.bio.indiana.edu/soft/molbio/readseq/java/readseq.jar
- 使用
# 提取數(shù)據(jù)未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
