软件：
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 #注释主程序，只能通过发邮件获取（试试http://www.openbioinformatics.org/annovar/download/0wgxR2rIVP/annovar.latest.tar.gz）
建库
mkdir ricedb
数据包括：
all.con #基因组序列
all.gff3 #MSU的v7.0版本组装的注释文件（上面两个文件http://rice.uga.edu/index.shtml）
gff文件会报错所以第一步要转换成gtf文件
gffread AsianRiice_MSU.gff3 -T -o AsianRice_MSU.gtf
gtf文件转换成GenePred文件，利用GtfToGenePred工具，这里注意“-genePredExt”这个参数一定要加上
gtfToGenePred -genePredExt AsianRice_MSU.gtf Os_refGene.txt
获得的GenePred文件：生成转录组信息文件
../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 表示使用基于基因的注释 一般是默认的#-dbtype refGene 表示使用"refGene"类型的数据库-out jpbc.anno 表示输出以jpbc.anno为前缀的结果文件
#Note:默认使用的为geneanno和refGene,这两个参数可以省略。#out:输出结果文件的前缀。#build:指定基因组构建版本。
假如vcf报错将vcf转化annovar
### 转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结果文件。
ricedb/数据库存放地址
5.对于多样本合并的vcf文件注释
用convert2annovar.pl对vcf4文件进行转换成annovar的输入文件格式时，对于有两个及以上ALT碱基类型，不同的参数会导致不同的结果
比如2780070位点有两个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
#转换时只有1型的ALT碱基保留下来，而2型的ALT碱基被过滤掉，其snp总行数与输入文件总行数相同
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
#在进行转换时ALT碱基的两种类型都会保留下来，并增加新增加一行，其总数比输入文件行数增多
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
在转换时如果想保留其他列的信息可以添加参数 --includeinfo
接下来是提取相应区间的变异信息:
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
该代码是一个shell脚本，它的作用是根据"pav-gene-pos.txt"文件中的位置信息对"sample.anno.variant_function"文件进行过滤，并将满足要求的行保存到文件中。
cat pav-gene-pos.txt | while read chr start end id：将文件"pav-gene-pos.txt"的内容逐行读取，并将每行的第一列数据赋值给变量"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"的值。过滤后的结果保存到tmp1.txt文件中。
cat tmp1.txt |while read line1：读取文件"tmp1.txt"的内容，逐行赋值给变量"line1"。
echo $line1 >> res/"${id}".csv：将变量"line1"的内容追加到名字为"${id}.csv"的文件中，该文件位于res目录下。
使用上述步骤来处理文件"pav-gene-pos.txt"中的每一行，直到所有行都处理完毕。
请注意，脚本中的注释行（以"##"开头的行）没有被执行，可以将它们删除。
