軟件：
1gffread
conda install -c bioconda gffread
gffread NIP-T2T.gff3 -T -o nip_t2t.gtf
2gffread #gff3 to gtf
gtfToGenePred? #gtf to genePred (建庫需要的文件)（下載： wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64.v369/gtfToGenePred）
chmod+x gtfToGenePred
3annovar? #注釋主程序堡牡，只能通過發(fā)郵件獲仁闾А（試試http://www.openbioinformatics.org/annovar/download/0wgxR2rIVP/annovar.latest.tar.gz）
建庫
mkdir ? ?ricedb
數(shù)據(jù)包括：
all.con #基因組序列
all.gff3 #MSU的v7.0版本組裝的注釋文件（上面兩個(gè)文件http://rice.uga.edu/index.shtml）
gff文件會(huì)報(bào)錯(cuò)所以第一步要轉(zhuǎn)換成gtf文件
gffread AsianRiice_MSU.gff3 -T -o AsianRice_MSU.gtf
gtf文件轉(zhuǎn)換成GenePred文件，利用GtfToGenePred工具晤柄，這里注意“-genePredExt”這個(gè)參數(shù)一定要加上
gtfToGenePred -genePredExt AsianRice_MSU.gtf Os_refGene.txt
獲得的GenePred文件：生成轉(zhuǎn)錄組信息文件
../retrieve_seq_from_fasta.pl --format refGene --seqfle all.fa Os_refGene.txt --out Os_refGeneMrna.fa
完成擦剑。得到Os_refGeneMrna.fa。 ?Os_refGene.txt
注釋：
perl table_annovar.pl ricedb/ --vcfinput --outfile fnal --buildver Os
注意#目前以上自建庫只能基于gene注釋#-geneanno 表示使用基于基因的注釋 一般是默認(rèn)的#-dbtype refGene 表示使用"refGene"類型的數(shù)據(jù)庫-out jpbc.anno 表示輸出以jpbc.anno為前綴的結(jié)果文件
#Note:默認(rèn)使用的為geneanno和refGene,這兩個(gè)參數(shù)可以省略芥颈。#out:輸出結(jié)果文件的前綴惠勒。#build:指定基因組構(gòu)建版本。
假如vcf報(bào)錯(cuò)將vcf轉(zhuǎn)化annovar
### 轉(zhuǎn)annovar~/tools/annovar/convert2annovar.pl-format vcf4 sample.vcf>sample.annovar
###注釋命令$>/public/home/fengting/tools/annovar/annotate_variation.pl-buildver gcf-geneanno ?sample.annovar ?-outfile ?sv.annovar ?ricedb/
sv.annovar結(jié)果文件爬坑。
ricedb/數(shù)據(jù)庫存放地址
5.對(duì)于多樣本合并的vcf文件注釋
用convert2annovar.pl對(duì)vcf4文件進(jìn)行轉(zhuǎn)換成annovar的輸入文件格式時(shí)纠屋，對(duì)于有兩個(gè)及以上ALT堿基類型，不同的參數(shù)會(huì)導(dǎo)致不同的結(jié)果
比如2780070位點(diǎn)有兩個(gè)ALT堿基
Chr1 2472436 . G T 22.72 PASS AC=2;AF=0.333;AN=6;BaseQRankSum=-0.262;DP=22;Dels=0.00;ExcessHet=0.4576;FS=0.000;HaplotypeScore=1.8128;MLEAC=2;MLEAF=0.333;MQ=23.58;MQ0=4;MQRankSum=-0.954;QD=11.36;ReadPosRankSum=-0.244;SOR=2.179 GT:AD:DP:GQ:PL 0/0:15,1:16:18:0,18,242 0/0:4,0:4:3:0,3,40 1/1:0,2:2:6:55,6,0
Chr1 2780070 . A C,T 1764.90 PASS AC=4,2;AF=0.667,0.333;AN=6;DP=46;Dels=0.00;ExcessHet=3.0103;FS=0.000;HaplotypeScore=0.0000;MLEAC=4,2;MLEAF=0.667,0.333;MQ=60.00;MQ0=0;QD=32.94;SOR=0.874 GT:AD:DP:GQ:PL 1/1:0,21,0:21:63:885,63,0,885,63,885 1/2:0,10,1:11:7:397,37,7,360,0,357 1/2:0,10,4:14:99:508,148,118,360,0,348
Chr1 5693321 . T G 616.12 PASS AC=4;AF=0.667;AN=6;BaseQRankSu
類型1 --format vcf4old
(base) fyx@huangji-5885H-V5:~/biosoft/annovar$perl convert2annovar.pl bo_r_s_3_snp.vcf -format vcf4old
Chr1 2472436 2472436 G T hom 22.72 22 23.58 11.36
Chr1 2780070 2780070 A C hom 1764.90 46 60.00 32.94
Chr1 5693321 5693321 T G hom 616.12 18 60.00 28.57
#轉(zhuǎn)換時(shí)只有1型的ALT堿基保留下來盾计，而2型的ALT堿基被過濾掉售担，其snp總行數(shù)與輸入文件總行數(shù)相同
NOTICE: Read 730 lines and wrote 681 different variants at 681 genomic positions (681 SNPs and 0 indels)
NOTICE: Among 681 different variants at 681 positions, 0 are heterozygotes, 681 are homozygotes
NOTICE: Among 681 SNPs, 470 are transitions, 211 are transversions (ratio=2.23)
類型2 -format vcf4 -allsample -withfreq 等同于-format vcf4old --allallele
(base) fyx@huangji-5885H-V5:~/biosoft/annovar$perl convert2annovar.pl bo_r_s_3_snp.vcf -format vcf4 -allsample -withfreq
Chr1 2472436 2472436 G T 0.3333 22.72 2
Chr1 2780070 2780070 A C 0.6667 1764.90 14
Chr1 2780070 2780070 A T 0.6667 1764.90 14
Chr1 5693321 5693321 T G 0.6667 616.12 1
#在進(jìn)行轉(zhuǎn)換時(shí)ALT堿基的兩種類型都會(huì)保留下來，并增加新增加一行署辉，其總數(shù)比輸入文件行數(shù)增多
NOTICE: Finished reading 730 lines from VCF file
NOTICE: A total of 681 locus in VCF file passed QC threshold, representing 710 SNPs (481 transitions and 229 transversions) and 0 indels/substitutions
NOTICE: Finished writing allele frequencies based on 2130 SNP genotypes (1443 transitions and 687 transversions) and 0 indels/substitutions for 3 samples
在轉(zhuǎn)換時(shí)如果想保留其他列的信息可以添加參數(shù) --includeinfo
接下來是提取相應(yīng)區(qū)間的變異信息:
cat pav-gene-pos.txt | while read chr start end id
do
##awk '{if($3 == "'"$chr"'" && $3 == "gene" && $4 >= "'"${start}"'" && $5 <= "'"${end}"'")print } ' sample.anno.variant_function > tmp1.txt
awk '{if($3 == "'"$chr"'" && $4 >= "'"${start}"'" && $5 <= "'"${end}"'")print } ' sample.anno.variant_function > tmp1.txt
cat tmp1.txt |while read line1
do
? ? echo $line1 >> res/"${id}".csv
done
##awk '$1 == "1" && $3 == "gene" && $4 >= 10000 && $5 <= 200000 {print $0} ' Arabidopsis_thaliana.TAIR10.31.gff3 > out.gff
done
該代碼是一個(gè)shell腳本族铆，它的作用是根據(jù)"pav-gene-pos.txt"文件中的位置信息對(duì)"sample.anno.variant_function"文件進(jìn)行過濾，并將滿足要求的行保存到文件中哭尝。
cat pav-gene-pos.txt | while read chr start end id：將文件"pav-gene-pos.txt"的內(nèi)容逐行讀取哥攘，并將每行的第一列數(shù)據(jù)賦值給變量"chr"，第二列賦值給變量"start"材鹦，第三列賦值給變量"end"逝淹，第四列賦值給變量"id"。
awk '{if($3 == "'"$chr"'" && $4 >= "'"${start}"'" && $5 <= "'"${end}"'")print } ' sample.anno.variant_function > tmp1.txt：使用awk命令過濾文件"sample.anno.variant_function"中的行桶唐。條件是第三列等于變量"chr"的值创橄，并且第四列大于或等于變量"start"的值，第五列小于或等于變量"end"的值莽红。過濾后的結(jié)果保存到tmp1.txt文件中。
cat tmp1.txt |while read line1：讀取文件"tmp1.txt"的內(nèi)容，逐行賦值給變量"line1"安吁。
echo $line1 >> res/"${id}".csv：將變量"line1"的內(nèi)容追加到名字為"${id}.csv"的文件中醉蚁，該文件位于res目錄下。
使用上述步驟來處理文件"pav-gene-pos.txt"中的每一行鬼店，直到所有行都處理完畢网棍。
請(qǐng)注意，腳本中的注釋行（以"##"開頭的行）沒有被執(zhí)行妇智，可以將它們刪除滥玷。
