先安装依赖的包（要先装BioConductor）：
假设fasta文件名为: aligned_fasta.fasta 读取fasta文件，转化：
结果文件为out.phy
注意:生成out.phy里,第一列序列名和第二列序列只有一个空格,而mcmctree要求两个以上,所以需要人工加多一个空格.
也可以用shell命令生成:
目录
vcfR 可以直接读取vcf格式的数据。如果同时读取参照序列fasta格式的序列文件和gff格式文件的注释文件还可以获取更完整的信息(此步骤并非必须，可以只读取vcf数据)。在此处便于重复用到了 pinfsc50 包。这个包里是植物致病微生物的基因序列测序结果。包含了一个vcf文件，一个fasta文件和一个gff文件。
这里用到参照序列的数据。
当这些数据被读取到内存的时候就可以开始对染色体名字或者其它一些东西进行修改了。由于 vcfR 更擅长对的单独染色体进行分析，所以当你的基因过大或者有很多样本的时候，建议对数据进行拆分。
读取完数据以后就可以建立 chromR ，来对数据进行详细的分析。
首先对数据进行初步的可视化,
我们在上面的图里得到很多信息，比方说测序深度(DP)的峰在500，但是拖着尾巴，这个尾巴表示数据里包含着CNV信息。然后比对质量(MQ)的峰值在60,于是我们可以以60为中心对数据进行过滤。
使用 masker 可以对数据进行过滤标记。然后可视化过滤以后的数据。
是不是顺眼多了。当然我们也可以看一下SNP的分布情况。注意右下角的图。
用 chromoqc() 可以对数据进行更完整的可视化。包括外显子内含子的分布，GC含量的分布等等。
最后可以用函数 write.vcf() 把数据输出成新的vcf文件。
计算β多样性指数需要用到phyloseq包。它的安装方式不同于简单的install.packages（“phyloseq”）
有两种方法可以安装
1.先安装BiocManager
install.packages("BiocManager")
library("BiocManager")
BiocManager::install("phyloseq")
library("phyloseq")
2.source("https://bioconductor.org/biocLite.R")
biocLite("phyloseq")
#安装phyloseq
library("phyloseq")
安装并加载了phyloseq包后，开始读取数据，前面计算α多样性，用到的是read.table……
qiimedata <- import_qiime(otufilename = "feature-table.taxonomy.txt", mapfilename = "mapping_file.txt", treefilename = "tree.rooted.nwk", refseqfilename = "dna-sequences.fasta")
#读取数据，参数都是文件名，注意加后缀
#otufilename指定out表格，mapfilename指定map文件（分组数据)
#treefilename指定有根进化树文件
#refseqfilename指定代表序列文件
otu<-qiimedata@[email protected]
#从qiimedata里面提取otu
sum_of_otus<-colSums(t(otu))
#t_转置,colsums计算列的和,即计算各个otu检测到的总序列数，为了筛掉一些总序列数过低的otu（可能是测序错误）
sum_of_otus
#查看otu总序列数
selected_otu<-names(sum_of_otus)[sum_of_otus>10]
#获取总序列数大于10的otu id
sub_qiimedata <- prune_taxa(selected_otu, qiimedata)
#筛选总序列数大于10的otu的phyloseq数据
weighted_unifrac<-distance(sub_qiimedata,method = 'wunifrac')
#计算样本间加权unifrac
unweighted_unifrac<-distance(sub_qiimedata,method = 'unifrac')
#计算样本间非加权unifrac
bray_curtis <- distance(sub_qiimedata, method='bray')
write.table(as.matrix(bray_curtis),"bray_curtis.txt",sep = '\t',quote = FALSE,col.names = NA)
#保存距离矩阵
#计算样本间Bray-Curtis距离矩阵，method 可选" wunifrac ", " unifrac " ，"jaccard"等
pcoa_of_bray_curtis<-ordinate(physeq=sub_qiimedata,distance = 'bray',method = "PCoA")
#基于Bray-Curtis距离矩阵的PCoA排序分析
p<-plot_ordination(sub_qiimedata, pcoa_of_bray_curtis, type="samples", color="Group1",shape = "Group1")
#将PCoA排序分析结果可视化
library("ggplot2")
p<-p+ scale_colour_manual(values=c("#DC143C","#808000","#00CED1")) + geom_point(size=2) +ggtitle("PCoA of Bray-Curtis distance")+theme(text = element_text(size = 15))
#修改图形大小,ggtitle加标题,stat_ellipse加椭圆
#用scale_colour_manual(values=c())自定义颜色，可查颜色的16进制对照表
p
nmds_of_bray_curtis<-ordinate(physeq=sub_qiimedata,distance = 'bray',method = "NMDS")
#基于Bray-Curtis距离矩阵的NMDS排序分析
p1<-plot_ordination(qiimedata, nmds_of_bray_curtis, type="samples", color="Group1")
#将NMDS排序分析结果可视化
# color=“Group1”指定不同分组的点染不同颜色
p1
p1<-p1+ geom_point(size=3) +ggtitle("NMDS of Bray-Curtis distance") + stat_ellipse()+theme(text = element_text(size = 15))
#对图片进行适当修饰， stat_ellipse()加椭圆， ggtitle()加标题
ggsave(plot = p1,“nmds_of_bary_curtis.pdf",dpi = 300,width
PCoA中的两个点距离，接近β多样性指数
PCA(Principal Components Analysis)即主成分分析，也称主分量分析或主成分回归分析法，首先利用线性变换，将数据变换到一个新的坐标系统中然后再利用降维的思想，使得任何数据投影的第一大方差在第一个坐标(称为第一主成分)上，第二大方差在第二个坐标(第二主成分)上。这种降维的思想首先减少数据集的维数，同时还保持数据集的对方差贡献最大的特征，最终使数据直观呈现在二维坐标系。
PCoA(Principal Co-ordinates Analysis)分析即主坐标分析，可呈现研究数据相似性或差异性的可视化坐标，是一种非约束性的数据降维分析方法，可用来研究样本群落组成的相似性或相异性。它与PCA类似，通过一系列的特征值和特征向量进行排序后，选择主要排在前几位的特征值，找到距离矩阵中最主要的坐标，结果是数据矩阵的一个旋转，它没有改变样本点之间的相互位置关系，只是改变了坐标系统。两者的区别为PCA是基于样本的相似系数矩阵(如欧式距离)来寻找主成分，而PCoA是基于距离矩阵(欧式距离以外的其他距离)来寻找主坐标。
NMDS图中两个点的距离的排序，接近β多样性指数的排序
