生物信息学常用软件

持续输出中。。。。。。

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

有一下这些文件：

1
2
3
4
5
6
7
8
9
10
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（所有的变异）文件是这样的：

1
2
3
4
5
6
7
8
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（外显子区域的变异）文件是这样的：

1
2
3
4
5
6
7
8
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 Analyses

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分析：

1
2
3
4
5
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 Analysis

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值最小。

1
2
3
4
5
6
7
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
1
2
3
4
5
6
7
8
9
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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
###
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 下载安装

1
2
3
4
5
6
7
8
9
10
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 使用案例

1
2
3
4
5
6
7
8
9
10
# 启动
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日
许可协议