生物信息学常用软件

持续输出中。。。。。。

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邮箱申请下载。下载完成后直接解压即可使用。

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文件：

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

然后提取基因的mRNA序列：

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文件进行转换：

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后面跟的是前两部使用到的前缀：

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

输出文件有这些：

-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安装

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

2.4 构建索引

bwa index NDH108.fa

2.5 序列比对

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进行安装：

mamba install bioconda::gatk4

3.4 构建索引

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

生成的文件是NDH108.dict.

3.5 变异检测

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，运行完成后会输出距离矩阵和系统发育树。

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文件时，输入文件需要是压缩过的。

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文件

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 使用案例

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

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抽样

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
  持续更新中。。。。。。