變異注釋
annovar 和 SnpEff 都是比較常用的 SNP 和 INDEL 注釋軟件嘱支。
annovar
ANNOVAR是一個(gè)高效的Perl 命令注釋工具命咐，可實(shí)現(xiàn)對(duì)來(lái)自下一代測(cè)序 （NGS） 數(shù)據(jù)的遺傳變異進(jìn)行高通量功能注釋和過(guò)濾沸毁，支持包括VCF在內(nèi)的多種輸入和輸出文件格式焚挠。
主要包含三種不同的注釋方法：gene-based, region-based 和filter-based制肮》硇基于基因的注釋（Gene-based Annotation）揭示variant與已知基因直接的關(guān)系以及對(duì)其產(chǎn)生的功能性影響，基于區(qū)域的注釋（Region-based Annotation）揭示variant 與不同基因組特定段的關(guān)系烦租，例如：它是否落在轉(zhuǎn)錄因子結(jié)合區(qū)域等延赌，基于篩選的注釋（Filter-based Annotation）則分析變異位點(diǎn)是否位于指定的數(shù)據(jù)庫(kù)中，比如dbSNP, 1000G叉橱，ESP 6500等數(shù)據(jù)庫(kù)挫以，計(jì)算SIFT/PolyPhen/LRT/MutationTaster/MutationAssessor等
- 數(shù)據(jù)：基因組文件（genome）、基因組注釋文件（gtf或gff）窃祝、VCF文件屡贺、INDEL文件
- 軟件：annovar,gtfToGenePred（gtf 文件做輸入）,gff3ToGenePred(gff3 文件做輸入)
# step1 轉(zhuǎn)gtf或gff成refGene格式并提取轉(zhuǎn)錄本的序列
genome.gtf
genome.fasta
## 轉(zhuǎn)gtf成refGene格式
gtfToGenePred -genePredExt genome.gtf genome_refGene.txt
##如果提供的為gff文件，使用下面命令
gff3ToGenePred genome.gff genome_refGene.txt
## 提取序列文件
/home/software/annovar/retrieve_seq_from_fasta.pl --format refGene --seqfile genome.fasta genome_refGene.txt --out genome_refGeneMrna.fa
# step2 annovar注釋
## 轉(zhuǎn)vcf成表格格式（軟件需求）
/home/software/annovar/convert2annovar.pl \
-format vcf4 \ #輸入文件格式
-allsample \ #allsample模式
-withfreq \ #輸出alt頻率
genome.vcf \ #輸入文件
> genome.annovar.input #輸出文件
## 進(jìn)行變異注釋 如果需要所有信息锌杀，添加-separate
/home/software/annovar/annotate_variation.pl \
-geneanno \ #基于基因進(jìn)行注釋
--neargene 2000 \ #確定gene上下游范圍
-buildver genome \ #基因組數(shù)據(jù)庫(kù)名稱
-dbtype refGene \ #數(shù)據(jù)庫(kù)類型
-outfile all.anno \ #輸出文件前綴
-exonsort \ # 結(jié)果按exon進(jìn)行排序
genome.annovar.input \ # 輸入文件
./ # 輸入數(shù)據(jù)庫(kù)路徑
# 輸出文件為：
## all.anno.exonic_variant_function
## all.anno.variant_function
# 基本統(tǒng)計(jì)
cut -f 1 all.anno.variant_function |sort | uniq -c > all.anno.variant_function.stat
cut -f 2 all.anno.exonic_variant_function |sort |uniq -c > all.anno.exonic_variant_function.stat
variant_function注釋結(jié)果甩栈，優(yōu)先級(jí)排序?yàn)?exonic = splicing > ncRNA > UTR5/UTR3 > intron > upstream/downstream > intergenic
INDEL注釋直接將腳本內(nèi)相應(yīng)snp文件替換成indel文件即可
SnpEff
SnpEff 軟件使用時(shí)需要進(jìn)行基因組數(shù)據(jù)庫(kù)構(gòu)建，這個(gè)軟件已經(jīng)預(yù)先構(gòu)建好了超過(guò) 2500 個(gè)物種的數(shù)據(jù)庫(kù), 如果使用的基因組版本在預(yù)先構(gòu)建的數(shù)據(jù)庫(kù)里存在糕再，可以直接下載使用量没。
如果研究的物種沒有預(yù)先建好的數(shù)據(jù)庫(kù)，可以自己手動(dòng)來(lái)構(gòu)建突想。
# 注釋
java -jar snpEff.jar -c ./snpEff.config \ #指定數(shù)據(jù)庫(kù)
-ud 2000 \ #指定上下游范圍
-csvStats test.csv \ #生成表格的統(tǒng)計(jì)文件
-htmlStats test.html \ #生成網(wǎng)頁(yè)的統(tǒng)計(jì)文件
-o vcf \ #輸出文件vcf格式
ric \ #數(shù)據(jù)庫(kù)名稱
genome.vcf > genome.ann.vcf
## 輸出文件文件 genome.ann.vcf
snpEff 變 異 類 型 說(shuō) 明 https://pcingola.github.io/SnpEff/se_inputoutput/#effect-prediction-details
對(duì)區(qū)域進(jìn)行注釋
BSA 分析得到的結(jié)果是基因組區(qū)域殴蹄，需要查看區(qū)域內(nèi)有哪些基因，可以使用 bedtools 進(jìn)行查看猾担。
# 準(zhǔn)備gff或gtf文件袭灯，以及需要查看的區(qū)域bed文件
genome.gtf
genome_region.bed
# 進(jìn)行區(qū)域內(nèi)基因查找
bedtools intersect -wo -a genome_region.bed -b genome.gtf | awk '$6=="gene"' > genome_region.anno
