👨‍👩‍👧‍👦LorMe大合成菌群高通量测序及生物信息学分析指南

LorMe大合成菌群由来自100个不同种的植物有益菌组成。完成试验并提取DNA后，根据试验需求进行相应送测。

LorMe大合成菌群高通量测序及生物信息学分析指南可以通过钉钉以及百度网盘下载获取：

钉钉网页链接：/全员文件夹/LorMe统计/LorMe大合成菌群高通量测序及生物信息学分析指南/

百度网盘链接：https://pan.baidu.com/s/1V-bIpTFtzr1koJ1X3SLC2w?pwd=1234 提取码: 1234

送测指南

纯培养推荐使用低深度的宏基因组测序，成本与二代扩增子测序差不多；土壤或根际样品依据成本与试验需求（详情请见“常见问题”模块）选择二代（便宜）、三代（略贵）或宏基因组测序（较贵）。
送测请联系：凌恩生物 程经理
1.扩增子测序
引物：二代选择515F_806R,全长选择27F_1492R
测序深度：通常30000，如果有特殊需求可增加测序数据量
指定流程：LorMe合成菌群扩增子测序标准流程
数据处理（丰度计算方法）：ASV
2. 宏基因组测序
请根据实际情况咨询程经理并与指导老师和汪宁祺进行确认。

测序结果文件校对

1）保存原始数据。获得测序结果后，首先妥善保存原始文件的压缩包。

2）校验测试质量。进入00.data文件夹，内含测序数据的统计文件，分别查看并校验①'cleanData.stat'中,所有样品的sequence分布及最低sequence数是否符合送测时的'测序深度'。②unclassify率。合成菌群试验理论上所有物种都应该是合成菌群内物种，根据送测内容查看对应的未注释到物种比例。通常，纯培养未注释率应在2%以下；盆栽土壤由于实验本身具有未封闭性或浇水时无法避免杂菌污染，未注释率会略高，具体请联系汪宁祺进行核验。如出现未注释率异常的情况，请检查样品是否污染/取样是否规范，如无以上情况，再联系公司进行流程比对参数的微调以降低未注释率。

3）提取对应文件进行本地处理。在测序质量符合规范的情况下，可进行选择ASV表进行本地处理。进入01.ASV文件，选择去除unclassify后再抽平的'asv_table_tax'表。

本地R语言处理

进入preprocessing文件夹进行数据预处理

1）数据编辑。打开'data'文件夹，参照'group_example.txt'文件编辑你的分组文件，并将其与asv表一起放入'data'文件夹中。

本地文件说明

-----preprocessing.Rproj  #预处理R项目文件
-----example_script.R     #预处理示例脚本
-----LorMe_0.6.0.tar.gz   #LorMe 0.6.0版本安装包

--Tree文件夹(可直接取用进化树及可视化文件，也可用于iTOL再加工)
-----strain.tre                        #菌株水平进化树 
-----species.tre                       #种水平进化树
-----iTOL_treecolors-phylum_strain.txt #菌株水平iTOL门映射注释
-----iTOL_treecolors-phylum-species.txt#种水平iTOL门映射注释
-----iTOL_tree_strain.pdf              #菌株水平系统发育进化树
-----iTOL_tree_species.pdf             #种水平系统发育进化树
-----iTOL_labels-strain.txt            #菌株水平iTOL标签注释
-----iTOL_labels-species.txt           #种水平iTOL标签注释

--data文件夹
-----syncom.fasta          #合成菌群种水平（CSCID）序列
-----Syncom_summary.R      #合成菌群预处理函数源代码
-----tax_infor_species.txt #species水平物种信息表
-----tax_infor_strain.txt  #strain水平物种信息表
-----group_example.txt     #处理的分组信息示例文件
-----asv_table_example.txt #合成菌群测序获得的asv示例文件
-----autoplot.R            #用于数据拓展的可视化APP源代码

--Figure文件夹(数据拓展跑出来的示例图片)

2）预处理。在本地打开“preprocessing.Rproj”文件，随后打开“example_script.R”脚本运行代码。

##功能与合成菌群信息加载####
source("./data/Syncom_summary.R") #加载预处理封装函数
tax_info=read.table("./data/tax_infor_species.txt",head=T,sep='\t') #读取species水平物种信息表
tax_info_strain=read.table("./data/tax_infor_strain.txt",head=T,sep='\t') #读取strain水平物种信息表

###合成菌群ASV表与分组处理信息加载####
asv=read.table("./data/asv_table_example.txt",head=T,sep="\t")  #读取合成菌群测序获得的asv示例文件
group_info=read.table("./data/group_example.txt",head=T,sep='\t')  #处理的分组信息示例文件

##数据处理####
Syncom_object=Syncom_summary(asvtable = asv,grouptable = group_info,syncom_tax = tax_info,strain_tax=tax_info_strain,outputstyle = "complete") #构建syncom对象

#数据提取示例
species=Syncom_object$Species              #种水平绝对丰度表
species_pct=Syncom_object$Species_percent  #种水平相对丰度表
species_taxonomy=Syncom_object$Species_taxonomy  #种水平完整注释

执行完成后，‘Syncom_object’即为预处理完的合成菌群对象。内含ASV及各分类水平绝对丰度、相对丰度和分类学注释的表，既可导出单独使用，也兼容LorMe包0.6.0以上版本进行数据拓展

3）数据拓展。预处理完成的‘Syncom_object’可兼容LorMe包0.6.0以上版本进行数据拓展(需本地安装文件夹中的LorMe0.6.0安装包；可继续运行代码，此处为中文注释讲解)

 #data expanding 数据拓展，一般用species水平分析
##package loading####  
source("./data/autoplot.R") #APP加载
pkg = c("shiny","shinydashboard","magrittr","shinyhelper","ggbeeswarm",'ggplot2',"ggpubr","agricolae","DescTools","multcompView","HH","coin") #依赖包加载
for(p in pkg){
  if(!requireNamespace(p, quietly = TRUE))
    install.packages(p)#, repos=site)
}
lapply(pkg, library, character.only = TRUE)

library(LorMe) #version 0.6.0+

community_alpha=Alpha_diversity_calculator(taxobj = Syncom_object,taxlevel = "Species") #计算alpha多样性

plot_maker() #实时可视化APP,调试完毕后请复制代码并关闭
#以下部分为由plot_maker() 自动生成的代码
my_col<-c('#35A585','#EAE48E','#006FB0','#CC78A6') 
names(my_col)<-c('T1','T2','T3','T4')#分组颜色设置
inputdata<- subset(community_alpha,Indexname %in% c("Shannon","ACE"))
facet_compare<-data.frame(compare=0,Letters=0,type=0,Mean=0,std=0,letterp=0,facetlabel=0)[0,]
for(i in unique(inputdata[,5])){
  sub_facet<-inputdata[which((inputdata[,5])==i),] 
  results<-auto_signif_test(data =sub_facet,treatment_col =2,value_col =6,prior = T)
  facet_compare<-rbind(facet_compare,data.frame(results$comparison_letters[,1:5],letterp=1.1*(results$comparison_letters %>% .[,'Mean'])+1.3*(results$comparison_letters %>% .[,'std']),facetlabel=i))
} #分面进行多重比较字母标记

colnames(facet_compare)[7]=colnames(inputdata)[5]
ggplot(inputdata,aes(x=as.factor(inputdata[,2]),y=inputdata[,6]))+
  scale_y_continuous(expand = c(0.05,0.01))+
  scale_color_manual(values=my_col)+
  scale_fill_manual(values=my_col)+
  theme_zg()+
  facet_wrap(~get(colnames(inputdata)[5]),scales ='free_y',strip.position ='top',as.table =T)+
  geom_boxplot(aes(fill=factor(inputdata[,2])),alpha=0.8,width=0.5,outlier.color =NA,color='#000000',linewidth=0.2)+
  labs(x='',fill='',y="")+
  geom_text(data=facet_compare,aes(x=as.factor(compare),y=letterp,label=Letters),size=6)+
  theme(legend.position = 'right') #可视化


community_composition=community_plot(taxobj = Syncom_object,taxlevel = "Species",n = 10,treatlocation = 2,replocation = 4,nrow = 1,treatorder = c("T1","T2","T3","T4")) #群落组成分析
community_composition$barplot            #群落组成柱状图
community_composition$areaplot           #群落组成area图
community_composition$alluvialplot       #群落组成冲积图
community_structure=structure_plot(taxobj = Syncom_object,taxlevel = "Species",treatlocation = 2,ellipse.level = .75) #群落结构分析
community_structure$PCA_Plot+
  scale_color_manual(values=my_col)    #PCA降维群落结构
community_structure$Pcoa_Plot+
  scale_color_manual(values=my_col)    #PCOA降维群落结构
community_structure$NMDS_Plot+
  scale_color_manual(values=my_col)    #NMDS降维群落结构

 #网络分析###
network_results=network_analysis(input = Syncom_object$ASV_percent,inputtype = 1,n = 10,threshold = .6,method = "spearman",display = T)  #网络分析
network_details=as.data.frame(network_results[1]) #节点信息表
network_adj=as.data.frame(network_results[2])     #邻接表
#以上两个表可无缝兼容Geghi与cytoscape进行可视化
##LorMe包其余功能可自由探索####

方法与材料写作模板

我们在文件夹中提供了一份供参考的'合成菌群数据库构建材料与方法英文模板'，请大家以此为参考进行写作，请勿原封不动照抄！

常见问题

1. 合成菌群物种个数问题。合成菌群菌液实际由来自100个非冗余物种的122个不同的植物有益菌株组成，在写作中，只需“合成菌群由100个物种构成”相关的描述即可

2. 数据分析选择的分类水平。按照每个species选择一株菌的理论模式，合成菌群数据分析通常直接基于“species”水平进行所有的数据分析，这样找到的关键物种直接对应明确的species，有助于进行后续的实验或其他组学验证。

3. 找到的关键species对应多个strain。使用合成菌群菌液时，共有来自100个species的122个strain。如需要验证的关键物种同时对应多个strain，则所对应的几个strain需同时进行验证，再根据其他生物学证据确定发挥作用的strain。具体请咨询相应指导老师和汪宁祺。

4. 以515F_806R为引物的二代扩增子测序分辨率问题。以515F_806R为引物二代扩增子测序选择V4区域进行扩增后比对，由于部分同属不同种的strain存在一模一样的V4区序列，导致二代扩增子测序无法准确分辨这些物种。如在分析完成后发现的关键species对应的strain存在与其他strain一致的V4区序列，则所有这些V4区序列一致的strain都需进行验证并根据其他生物学证据确定发挥作用的strain。具体请咨询相应指导老师和汪宁祺。全长测序和宏基因组测序则可有效解决以上问题（就是测序成本会比较高）。

5. 包含其他物种的合成菌群试验高通量测序。试验包含其他物种（如青枯菌、链霉菌、假单胞菌）时，需向测序公司提供相应物种的16s序列用以建库。首推该物种全基因组测序中抽取的16s序列，如没有，则使用全长引物进行pcr扩增获得该物种16s全长序列（全长引物27F_1492R在首尾部分会有几十个碱基的误差，准确度不如全基因组中抽取的16s序列）。送测时，请在样品单中说明需加入相应物种的序列用以建库。本地处理方面，需在'tax_infor_species.txt'和'tax_infor_strain.txt'两个文件中补充相应物种的信息，具体方式请咨询汪宁祺。

6. 合成菌群物种分类问题。 1). 当前我们以合成菌群全基因组测序合成菌群ANI的注释结果为准，确定了共100个非冗余的种（species）。如单独将菌株以全长引物扩增获得16s序列，再进行NCBI比对，可能会出现鉴定出来结果的第一位不一定是我们标注的物种，这种情况的出现原因同第四条（全长引物本身有几十个碱基的差异）。因此，使用我们提供的基于ANI抽取的16s序列构建的数据库以及在本地文件中提供的物种分类信息即可，我们已经过多方校验以保证每个物种的序列准确无误。 2). 由于NCBI对Tax的分类不断变化，部分物种的拉丁名分类注释已按照NCBI最新版本数据库更新。由于NCBI分类注释的数据库会保持更新，我们记录了每个菌株的TaxID（可在NCBI搜索获得相应TaxID的完整信息），以保证我们对各菌株的可追踪性。

技术支持

如有合成菌群各方面的疑惑或需要技术支持，请咨询：

实验部分——杨天杰

合成菌群全基因组和宏基因组应用、序列数据库构建与注释——张耀中

扩增子流程——杨欣润

数据校验、本地处理和数据拓展——汪宁祺