生物信息学常用软件

持续输出中。。。。。。

1 Annovar

1.1 参考文献

ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data

1.2 功能概述

对变异（SNP和INDLE等）进行注释。

1.3 下载安装

从https://www.openbioinformatics.org/annovar/annovar_download_form.php使用edu邮箱申请下载。下载完成后直接解压即可使用。

1
2
3

tar -zxf annovar.latest.tar.gz

cd annovar

有一下这些文件：

ll ~/software/annovar                                          
总计 532K
-rwxr-xr-x 1 lixiang lixiang 221K  6月  8  2020 annotate_variation.pl
-rwxr-xr-x 1 lixiang lixiang  42K  6月  8  2020 coding_change.pl
-rwxr-xr-x 1 lixiang lixiang 167K  6月  8  2020 convert2annovar.pl
drwxr-xr-x 2 lixiang lixiang 4.0K  6月  8  2020 example
drwxr-xr-x 3 lixiang lixiang 4.0K  6月  8  2020 humandb
-rwxr-xr-x 1 lixiang lixiang  19K  6月  8  2020 retrieve_seq_from_fasta.pl
-rwxr-xr-x 1 lixiang lixiang  42K  6月  8  2020 table_annovar.pl
-rwxr-xr-x 1 lixiang lixiang  21K  6月  8  2020 variants_reduction.pl

1.4 构建数据库

先使用gff3ToGenePred（下载完成后添加可执行权限即可：chmod +x gtfToGenePred）处理gff文件：

1	`gff3ToGenePred 01.data/genome/NDH108.gff3 05.annovar/NDH108_refGene.txt`

然后提取基因的mRNA序列：

1	`perl ~/software/annovar/retrieve_seq_from_fasta.pl --format refGene --seqfile 01.data/genome/NDH108.fa 05.annovar/NDH108_refGene.txt -outfile 05.annovar/NDH108_refGeneMrna.fa`

需要注意的是这两步得到的输出文件的前缀需要完全一致，就比如上面的NDH108，而且前缀的后面需要跟上下划线，refGene.txt和refGeneMrna.fa这两个名称必须这样写。

1.5 变异注释

首先对VCF文件进行转换：

1	`perl ~/software/annovar/convert2annovar.pl -format vcf4 -allsample -withfreq 04.vcf/G1.G2.merged.filtered.vcf > 04.vcf/G1.G2.merged.filtered.annovar.vcf`

然后进行注释，输入文件是上一部转换后的VCF文件，--buildver后面跟的是前两部使用到的前缀：

1	`perl ~/software/annovar/annotate_variation.pl --geneanno -dbtype refGene --buildver NDH108 04.vcf/G1.G2.merged.filtered.annovar.vcf 05.annovar -out 04.vcf/vcf.annov`

输出文件有这些：

1
2
3

-rw-rw-r-- 1 lixiang lixiang 2.1M  5月 29 11:59 vcf.annov.exonic_variant_function
-rw-rw-r-- 1 lixiang lixiang 1.1K  5月 29 11:59 vcf.annov.log
-rw-rw-r-- 1 lixiang lixiang 110M  5月 29 11:59 vcf.annov.variant_function

其中vcf.annov.variant_function（所有的变异）文件是这样的：

downstream      NDH01G00050(dist=428)   Chr01   86174   86174   -       TA      hom     626.6   15
intergenic      NDH01G00050(dist=1423),NDH01G00060(dist=18688)  Chr01   87169   87169   -       ATT     hom     775.61  19
intergenic      NDH01G00050(dist=2364),NDH01G00060(dist=17740)  Chr01   88110   88117   TATATATA        -       hom     857.3   21
intergenic      NDH01G00050(dist=2731),NDH01G00060(dist=17380)  Chr01   88477   88477   T       C       hom     1297.06 27
intergenic      NDH01G00050(dist=2908),NDH01G00060(dist=17203)  Chr01   88654   88654   G       A       hom     1844.06 41
intergenic      NDH01G00050(dist=3304),NDH01G00060(dist=16800)  Chr01   89050   89057   ATATATAT        -       het     462.02  12
intergenic      NDH01G00050(dist=3306),NDH01G00060(dist=16800)  Chr01   89052   89057   ATATAT  -       het     462.02  12
intergenic      NDH01G00050(dist=5367),NDH01G00060(dist=14744)  Chr01   91113   91113   -       GATGATGATGATGATGATGAT   hom     519.02  12

vcf.annov.exonic_variant_function（外显子区域的变异）文件是这样的：

line45  synonymous SNV  NDH01G00110:NDH01G00110.1:exon5:c.A2199G:p.R733R,       Chr01   174388  174388  A       G       hom     1406.06 36
line53  nonsynonymous SNV       NDH01G00140:NDH01G00140.1:exon1:c.G376A:p.V126I,NDH01G00140:NDH01G00140.2:exon1:c.G376A:p.V126I,        Chr01   206128  206128      G       A       hom     1673.06 44
line120 synonymous SNV  NDH01G00170:NDH01G00170.2:exon1:c.G33A:p.A11A,NDH01G00170:NDH01G00170.1:exon1:c.G33A:p.A11A,    Chr01   262077  262077  G       A           het     1459.64 87
line121 synonymous SNV  NDH01G00170:NDH01G00170.2:exon1:c.C153A:p.V51V,NDH01G00170:NDH01G00170.1:exon1:c.C153A:p.V51V,  Chr01   262197  262197  C       A           het     1379.64 79
line122 nonsynonymous SNV       NDH01G00170:NDH01G00170.2:exon1:c.C187G:p.R63G,NDH01G00170:NDH01G00170.1:exon1:c.C187G:p.R63G,  Chr01   262231  262231  C           G       het     1756.64 83
line123 nonsynonymous SNV       NDH01G00170:NDH01G00170.2:exon1:c.G191C:p.G64A,NDH01G00170:NDH01G00170.1:exon1:c.G191C:p.G64A,  Chr01   262235  262235  G           C       het     1756.64 83
line124 nonsynonymous SNV       NDH01G00170:NDH01G00170.2:exon1:c.C265A:p.R89S,NDH01G00170:NDH01G00170.1:exon1:c.C265A:p.R89S,  Chr01   262309  262309  C           A       het     1443.64 76
line125 synonymous SNV  NDH01G00170:NDH01G00170.2:exon1:c.T267G:p.R89R,NDH01G00170:NDH01G00170.1:exon1:c.T267G:p.R89R,  Chr01   262311  262311  T       G           het     1443.64 76

2 BWA

2.1 参考文献

Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM

2.2 功能概述

主要是序列比对。

2.3 下载安装

可以直接从GitHub编译安装，也可以用Conda安装

1 2	`git clone https://github.com/lh3/bwa.git cd bwa; make`

2.4 构建索引

1	`bwa index NDH108.fa`

2.5 序列比对

1	`bwa mem -t 70 -M -P 01.data/genome/NDH108.fa 01.data/raw/G1_CleanData_R1.fastq.gz 01.data/raw/G1_CleanData_R2.fastq.gz > 03.mapping/G1.sam`

3 GATK

3.1 参考文献

Genomics in the cloud: using Docker, GATK, and WDL in Terra

3.2 功能概述

主要是对变异进行鉴定。

3.3 下载安装

建议使用Conda进行安装：

1	`mamba install bioconda::gatk4`

3.4 构建索引

1	`gatk CreateSequenceDictionary -R 01.data/genome/NDH108.fa`

生成的文件是NDH108.dict.

3.5 变异检测

1	`gatk --java-options "-Xmx60g" HaplotypeCaller -R 01.data/genome/NDH108.fa -I 03.mapping/P3_2.cleaned.fixed.group.DEDUP.bam -O 04.vcf/P3_2.vcf`

4 VCFDis

4.1 参考文献

VCF2Dis: an ultra-fast and efficient tool to calculate pairwise genetic distance and construct population phylogeny from VCF files

4.2 功能概述

使用VCF文件计算p-distance并构建群体的系统发育树。

4.3 下载安装

直接从GitHub上下载编译即可，也可以使用Docker等容器技术。

4.4 使用案例

测试使用的是Docker，运行完成后会输出距离矩阵和系统发育树。

1
2
3

docker run -it -v $(pwd):/data --name=vcf2dis registry.cn-shenzhen.aliyuncs.com/knight134/vcf2dis:v1.53e /bin/bash

VCF2Dis	-InPut	in.vcf.gz	-OutPut p_dis.mat

5 BCFtools

5.1 参考文献

Twelve years of SAMtools and BCFtools.

5.2 功能概述

处理VCF文件。

5.3 使用案例

5.3.1 合并VCF文件

在合并VCF文件时，输入文件需要是压缩过的。

1	`bcftools merge -o G1.G2.P3_1.P3_2.merged.vcf G1.vcf.gz G2.vcf.gz P3_1.vcf.gz P3_2.vcf.gz`

5.3.2 过滤VCF文件

1	`bcftools filter -i 'QUAL>30 && FORMAT/DP[*]>10' G1.G2.P3_1.P3_2.merged.vcf -o G1.G2.P3_1.P3_2.merged.filtered.vcf`

这个代码表示保留质量值大于30，而且所有样品的覆盖度都大于10的变异位点。

6 VCFtools

6.1 参考文献

The variant call format and VCFtools

6.2 功能概述

处理VCF文件。

6.3 使用案例

6.3.1 筛选SNP

1	`vcftools --vcf G1.G2.merged.filtered.vcf --remove-indels --recode --out G1.G2.merged.snp.vcf`

6.3.2 计算SNP密度

1	`vcftools --vcf G1.G2.merged.snp.recode.vcf --SNPdensity 10000 --out snp_density`

输出文件snp_density.snpden是SNP密度文件。

6.3.3 筛选INDEL

1	`vcftools --vcf G1.G2.merged.filtered.vcf --keep-only-indels --recode --out G1.G2.merged.indel`

6.3.4 计算INDEL密度

1	`vcftools --vcf G1.G2.merged.indel.recode.vcf --SNPdensity 10000 --out indel_density`

输出文件indel_density.snpden`是SNP密度文件。

7 plink

7.1 参考文献

PLINK: A Tool Set for Whole-Genome Association and Population-Based Linkage Analyses61352-4?_returnURL=https%3A%2F%2Flinkinghub.elsevier.com%2Fretrieve%2Fpii%2FS0002929707613524%3Fshowall%3Dtrue)

7.2 功能概述

全基因组数据分析软件。

7.3 使用案例

7.3.1 群体PCA

PCA分析前需要过滤一下SNP的缺失率（max-missing）、最小次等位基因频率（maf）和连锁SNP，使用plink完成。输入文件为硬过滤后的plink格式文件，见SNP calling文章的vcf文件转plink文件。
缺失率（max-missing）和次等位基因频率（maf）分别设置为0.2及0.01，也可以根据需要调整（注意vcftools和plink的max-missing数值设置是相反的）。
下面的代码是群体遗传过滤VCF文件用于PCA分析：

plink --bfile SNP --geno 0.2 --maf 0.01 --make-bed --out SNP.filter
plink --bfile SNP.filter --indep-pairwise 50 1 0.5 -out SNP.filter.pruning
plink --bfile SNP.filter --extract SNP.filter.pruning.prune.in --make-bed --out SNP.filter.pruned

plink --thread-num 2 --bfile SNP.filter.pruned --make-grm-bin --out SNP.filter.pruned.grm

8 GCTA

8.1 参考文献

GCTA: A Tool for Genome-wide Complex Trait Analysis00598-7?_returnURL=https%3A%2F%2Flinkinghub.elsevier.com%2Fretrieve%2Fpii%2FS0002929710005987%3Fshowall%3Dtrue)

8.2 功能概述

全基因组复杂形状分析软件。

8.3 使用案例

8.3.1 PCA

使用plink过滤后的VCF文件进行PCA分析。

1	`gcta-1.94.1 --thread-num 2 --grm SNP.filter.pruned.grm --pca 4 --out SNP.filter.pruned.pca`

9 admixture

9.1 参考文献

Fast model-based estimation of ancestry in unrelated individuals

9.2 功能概述

群体结构分析。

9.3 使用案例

使用PCA分析中过滤LD后的文件作为输入。根据实际情况，我们需要测试预设的群体数K从2到自己的总群体数。如果不确定总群体数，可以从2开始测试，直到CV error值最小。

ln -s PCA/SNP.filter.pruned.bed
ln -s PCA/SNP.filter.pruned.bim
ln -s PCA/SNP.filter.pruned.fam

plink --bfile SNP.filter.pruned --recode 12 --out SNP.filter.pruned.subset

admixture -j10 --cv SNP.filter.pruned.subset.ped 2

10 seqkit

10.1 参考文献

SeqKit2: A Swiss army knife for sequence and alignment processing

10.2 功能概述

对fastq或者是fasta格式的数据进行操作。

10.3 使用案例

10.3.1 reads抽样

1 2	`seqkit sample -s 100 -p 0.25 --two-pass all.R1.fastq -o all.R1.new.fastq seqkit sample -s 100 -p 0.25 --two-pass all.R2.fastq -o all.R2.new.fastq`

-s表示的是随机数种子，-p表示的是抽取百分之多少的reads.

10.3.2 提取序列

ID存放在文件中：seqkit grep -f 五种方法都有的基因ID.txt NDH108.pep.fa > target.pep.fa
单个ID进行提取：seqkit grep -p "NDH01G36010.1" NDH108.pep.fa
11 Augustus
11.1 参考文献
Stanke M, Diekhans M, Baertsch R, et al. Using native and syntenically mapped cDNA alignments to improve de novo gene finding[J]. Bioinformatics, 2008, 24(5): 637-644.
11.2 功能概述
主要用于基因组注释等。
11.3 使用案例
11.3.1 模型训练

新建一个物种

1	`new_species.pl --species=Rice_35minicore_NLR_0.1`

Will create parameters for a EUKARYOTIC species!
creating directory /share/org/YZWL/yzwl_lixg/miniforge3/envs/Augustus_v.3.5.0/config/species/Rice_35minicore_NLR_0.1/ ...
creating /share/org/YZWL/yzwl_lixg/miniforge3/envs/Augustus_v.3.5.0/config/species/Rice_35minicore_NLR_0.1/Rice_35minicore_NLR_0.1_parameters.cfg ...
creating /share/org/YZWL/yzwl_lixg/miniforge3/envs/Augustus_v.3.5.0/config/species/Rice_35minicore_NLR_0.1/Rice_35minicore_NLR_0.1_weightmatrix.txt ...
creating /share/org/YZWL/yzwl_lixg/miniforge3/envs/Augustus_v.3.5.0/config/species/Rice_35minicore_NLR_0.1/Rice_35minicore_NLR_0.1_metapars.cfg ...
The necessary files for training Rice_35minicore_NLR_0.1 have been created.
Now, either run etraining or optimize_parameters.pl with --species=Rice_35minicore_NLR_0.1.
etraining quickly estimates the parameters from a file with training genes.
optimize_augustus.pl alternates running etraining and augustus to find optimal metaparameters.

构建训练用的数据
使用的脚本是gff2gbSmallDNA.pl
拆分数据集
使用脚本把上一步准备的数据拆分为训练集和测试集，使用的脚本是randomSplit.pl
训练模型
使用etraining对模型进行训练
优化模型
使用optimize_augustus.pl优化模型

测试模型
把训练好的模型用来测试前面准备的测试集：

1	`augustus --species=Rice_35minicore_NLR_0.1 02.augustus/0.1_NLR_training_set.gb.test > 04.valid/0.1_prediction_result.gtf`

输出的文件的最后又预测结果的相关信息：

###
a-posteriori probability of viterbi path
----------------------------------------
a-posteriori probability of correct path

1008 times were the paths equally likely (identical).
sorted quotients of the rest:

0 quotients were between 1 and 10


*******      Evaluation of gene prediction     *******

---------------------------------------------\
                 | sensitivity | specificity |
---------------------------------------------|
nucleotide level |       0.988 |       0.911 |
---------------------------------------------/

----------------------------------------------------------------------------------------------------------\
           |  #pred |  #anno |      |    FP = false pos. |    FN = false neg. |             |             |
           | total/ | total/ |   TP |--------------------|--------------------| sensitivity | specificity |
           | unique | unique |      | part | ovlp | wrng | part | ovlp | wrng |             |             |
----------------------------------------------------------------------------------------------------------|
           |        |        |      |                866 |                715 |             |             |
exon level |   2893 |   2742 | 2027 | ------------------ | ------------------ |       0.739 |       0.701 |
           |   2893 |   2742 |      |  488 |   49 |  329 |  522 |   60 |  133 |             |             |
----------------------------------------------------------------------------------------------------------/

----------------------------------------------------------------------------\
transcript | #pred | #anno |   TP |   FP |   FN | sensitivity | specificity |
----------------------------------------------------------------------------|
gene level |  1139 |  1008 |  556 |  583 |  452 |       0.552 |       0.488 |
----------------------------------------------------------------------------/

------------------------------------------------------------------------\
            UTR | total pred | CDS bnd. corr. |   meanDiff | medianDiff |
------------------------------------------------------------------------|
            TSS |         40 |              0 |         -1 |         -1 |
            TTS |         66 |              0 |         -1 |         -1 |
------------------------------------------------------------------------|
            UTR | uniq. pred |    unique anno |      sens. |      spec. |
------------------------------------------------------------------------|
                |  true positive = 1 bound. exact, 1 bound. <= 20bp off |
 UTR exon level |          0 |              0 |       -nan |       -nan |
------------------------------------------------------------------------|
 UTR base level |          0 |              0 |       -nan |       -nan |
------------------------------------------------------------------------/
nucUTP= 0 nucUFP=0 nucUFPinside= 0 nucUFN=0
# total time: 246
# command line:
# /share/org/YZWL/yzwl_lixg/miniforge3/envs/Augustus_v.3.5.0/bin/augustus --species=Rice_35minicore_NLR_0.1 02.augustus/0.1_NLR_training_set.gb.test

12 edgeturbo

12.1 参考文献

12.2 功能概述

主要是用于下载国家生物信息中心的数据。

12.3 下载安装

mkdir ~/software/edgeturbo    
cd ~/software/edgeturbo

wget https://ngdc.cncb.ac.cn/ettrans/download/edgeturbo-client.linux.latest.cncb.tar.gz
tar -zxvf edgeturbo-client.linux.latest.cncb.tar.gz

# 编辑~/.zshrc文件
alias edge="/share/org/YZWL/yzwl_lixg/software/edgeturbo/edgeturbo-client/edgeturbo"

source ~/.zshrc

12.4 使用案例

# 启动
edge start

# 重启
edge restart

# 停止
egde stop

13 XtractPAV

13.1 参考文献

[]

13.2 功能概述

在多个基因组之间鉴定PAV.

13.3 使用案列

持续更新中。。。。。。

生物信息学

#生物信息学

生物信息学常用软件

https://lixiang117423.github.io/article/bioinftools/

作者

李详【Xiang LI】

发布于

2025年5月29日

许可协议

结构方程模型计算Fisher's C值上一篇

Claude给的基因组组装和注释的流程代码下一篇

生物信息学常用软件

1 Annovar

1.1 参考文献

1.2 功能概述

1.3 下载安装

1.4 构建数据库

1.5 变异注释

2 BWA

2.1 参考文献

2.2 功能概述

2.3 下载安装

2.4 构建索引

2.5 序列比对

3 GATK

3.1 参考文献

3.2 功能概述

3.3 下载安装

3.4 构建索引

3.5 变异检测

4 VCFDis

4.1 参考文献

4.2 功能概述

4.3 下载安装

4.4 使用案例

5 BCFtools

5.1 参考文献

5.2 功能概述

5.3 使用案例

5.3.1 合并VCF文件

5.3.2 过滤VCF文件

6 VCFtools

6.1 参考文献

6.2 功能概述

6.3 使用案例

6.3.1 筛选SNP

6.3.2 计算SNP密度

6.3.3 筛选INDEL

6.3.4 计算INDEL密度

7 plink

7.1 参考文献

7.2 功能概述

7.3 使用案例

7.3.1 群体PCA

8 GCTA

8.1 参考文献

8.2 功能概述

8.3 使用案例

8.3.1 PCA

9 admixture

9.1 参考文献

9.2 功能概述

9.3 使用案例

10 seqkit

10.1 参考文献

10.2 功能概述

10.3 使用案例

10.3.1 reads抽样

10.3.2 提取序列

11 Augustus

11.1 参考文献

11.2 功能概述

11.3 使用案例

11.3.1 模型训练

12 edgeturbo

12.1 参考文献

12.2 功能概述

12.3 下载安装

12.4 使用案例

13 XtractPAV

13.1 参考文献

13.2 功能概述

13.3 使用案列