[2021数学建模美赛A题] 如何获得真菌耐湿性数据、画出题目的图C和图A
首先列出可能用到的文章：
A题题目《2021_MCM_Problem_A》（以下简称题目A）
题目A出题参考《A trait-based understanding of wood decompositionby fungi》（https://www.pnas.org/content/pnas/117/21/11551.full.pdf）（以下简称论文A）
论文A的参考资料25《Consistent trade-offs in fungal trait expression across broad spatial scales》（以下简称论文25）
论文A的参考资料34《Diversity begets diversity in competition for space》（以下简称论文34）
然后列出能够找到的数据：
论文A的附录 (https://www.pnas.org/highwire/filestream/926186/field_highwire_adjunct_files/0/pnas.1909166117.sapp.pdf)（以下简称附录A）
附录A的参考资料9 (http://ncf.sobek.ufl.edu/AA00026436)
论文25 "Data availability" 章节提供的论文所用数据和代码的Github地址 (https://github.com/dsmaynard/fungal_biogeography)
论文34 "Data availability" 章节提供的论文所用数据和代码的Github地址 (https://github.com/dsmaynard/diversity_begets_diversity)
附录A提供了S1-S5五张表，其中S3和S4给出了34个菌种的木材分解率和菌丝延伸率数据
使用S3和S4的数据就可以画出题目的图C了（忽略不确定度）
图A的绘制就比较麻烦了，需要找到 moisture trade-off 数据。
题目A中对图A横轴的说明如下：
耐湿性（每种分离物的竞争排名及其水分位宽度的差异，均缩放至[0,1]）
moisture tolerance (difference of each isolate’s competitive ranking and their moisture niche width, both scaled to [0,1])
论文A中对图A横轴的说明如下：
水分权衡（A中的x轴）是每个分离物的竞争排名与其水分生态位宽度之间的差异，两者均按比例缩放至[0,1]，如Maynard等人(25)所计算。
The moisture trade-off (x axis in A) is the difference between each isolate’s competitive ranking and their moisture niche width, both scaled to [0,1], as calculated by Maynard et al. (25).
浏览一下论文25，可以看到 Fig. 4(a) 的纵轴为 "Dominance-tolerance" ，与图A的横轴 "moisture trade-off" 一样，取值范围均为 [-1,1] 。论文25对左图(a)的纵轴的说明为“耐热性权衡” (thermal dominance-tolerance trade-off) 。
Partial regression plot for PCA1 versus the thermal dominance-tolerance trade-off, with the solid line denoting the mean regression trend and the shaded region denoting the 95% confidence interval for the mean (n = 37 biologically independent isolates).
论文25提供了数据与代码的下载，其中 Fungal_climate_data.csv 提供了37个菌种的竞争排名 Competitive ranking(ranking)、温度生态位宽度 Temperature niche width(temp.niche.width)、水分生态位宽度 Moisture niche width(water.niche.width) 等数据，Fungal_trait_data.csv 除提供上述三种数据外，还有生长率 Growth rate(rate.0.5) 、菌丝密度 Hyphal density(density) 数据。需要注意，两张表格中的竞争排名 Competitive ranking(ranking) 数据并不一致。
由此推测，论文作者使用“温度生态位宽度”数据计算出“耐热性”，进而绘制出 Fig. 4(a) ，那么改用“水分生态位宽度”数据应该就能得到“耐湿性”。为了得知论文作者的计算过程，阅读论文25提供的R语言代码，可以找到 "# main temperature plot" 注释的部分和 av_gg 、scale01 等函数。
通过网络搜索对代码中一些R语言关键字的介绍、对R语言有了简单了解后，我安装了R和RStuido以运行代码。打开文件 Regression_and_maps.R ，从头执行至下述语句后，RStuido输出了论文25中的 Fig. 4 。
# main temperature plot
gg3<-plot_fig(mod=fit3.s,climate.pca=climate.pca, bstack.o=bstack.o, kc0=kc0, forest_mask=forest_mask, lake_fortify=lake_fortify, ocean_fortify=ocean_fortify,scale_func=scale11)
gg3av<-av_gg(fit3.s,"PC3.p",xlab="Climate PCA 1")
plot_grid(gg3av,gg3,rel_widths=c(.7,1),labels="AUTO")
重读代码，筛选出得到左图(A)所必须的代码行如下：
comb<-read_csv("../fungi_data/Fungal_climate_data.csv") %>% data.frame()
# scale the variables to be within [0,1]
comb$temp01<-scale01(comb$temp.niche.width)
comb$water01<-scale01(comb$water.niche.width)
comb$tw01<-scale01(comb$temp01*comb$water01)
# create the trade-off variables
comb <- comb %>% mutate(dom_tol=ranking-tw01) %>%
mutate(tdom_tol=ranking-temp01) %>%
mutate(wdom_tol=ranking-water01) %>%
mutate(water01=sqrt(water01)) %>%
mutate(temp01=sqrt(temp01)) %>%
mutate(ranking=sqrt(ranking))
# main regression models
summary(fit3.s<-lm(tdom_tol~ PC2.p + PC3.p + PC2.t + a.gal + p.fla + p.ruf,data=comb))
（中略）
gg3av<-av_gg(fit3.s,"PC3.p",xlab="Climate PCA 1")
这样，变量gg3av中就存有左图(A)。位于文件 misc_functions.R 的函数 av_gg 的代码如下：
av_gg<-function(mod,var="PC3.p",xlab="Climate PCA 3",ylab="Dominance-tolerance",col_pallette=plasma){
a<-data.frame(avPlots(fit3.s,terms=var))
names(a)<-c("x","y")
a <- a %>% mutate(y=((scale01(y)-0.5)*2))
gg2<-ggplot(a,aes(x=x,y=y,col=y))+geom_point()+
geom_smooth(method=lm,aes(color=..y..), size=1.5) +
scale_color_gradientn(colors=rev(col_pallette(100)),name="")+
xlab(xlab)+ylab(ylab)+ theme_bw()+
theme(legend.position = "none")
return(gg2)
}
可知变量a中保存有绘图用到的xy数据，图像保存在变量gg2中并作为函数返回值。利用函数 av_gg 的代码，将 Regression_and_maps.R 中 gg3av<-av_gg 这一行替换为如下三行：
a<-data.frame(avPlots(fit3.s,terms="PC3.p"))
names(a)<-c("x","y")
fit3 <- a %>% mutate(y=((scale01(y)-0.5)*2))
运行这三行代码，就得到了保存在变量fit3中的xy数据：
保存数据到csv文件中，使用Python绘图：
观察 Regression_and_maps.R 中如下几行：
comb$temp01<-scale01(comb$temp.niche.width)
comb$water01<-scale01(comb$water.niche.width)
（中略）
mutate(tdom_tol=ranking-temp01) %>%
mutate(wdom_tol=ranking-water01) %>%
（中略）
summary(fit3.s<-lm(tdom_tol~ PC2.p + PC3.p + PC2.t + a.gal + p.fla + p.ruf,data=comb))
summary(fit4.s<-lm(wdom_tol~ PC2.p + PC3.p + PC2.t + a.gal + m.tre + p.fla + p.ruf,data=comb))
可以发现，fit3.s与“耐热性”有关，fit4.s与“耐湿性”有关。把之前的三行代码修改如下并运行：
a<-data.frame(avPlots(fit4.s,terms="PC3.p"))
names(a)<-c("x","y")
fit4 <- a %>% mutate(y=((scale01(y)-0.5)*2))
这里的y列数据应该就是“耐湿度”(moisture tolerance)了。注意到论文25提供的是37个菌种的数据，而附录A中只有34种，经过对菌种名称的比对，需要去掉编号16、19、20的菌种数据。
观察图A，只有一种颜色的点。它的横轴数据应该是上面的y列数据，但它的纵轴数据（对数）“分解率”(decomposition rate) 不确定是哪个温度条件下的。把三个温度下的散点图都绘制出来看一下：
观察最右边两个点，似乎图A用的是10℃条件下的数据，但这两个点的纵坐标小于2，与图A不符。查看变量fit4，“耐湿度(moisture tolerance)”为最大值1的菌种编号36，名称为 "Tyromyces chioneus HHB11933 B10F"。图A中该点纵坐标约为2.6，e^2.6≈13.46。查看附录A的S3：
所以图A的y轴数据应该是附录A的S3的最右边一列“几何平均值(geometric mean)”。再次绘图：
还是不对。再次阅读代码：
a<-data.frame(avPlots(fit4.s,terms="PC3.p"))
names(a)<-c("x","y")
fit4 <- a %>% mutate(y=((scale01(y)-0.5)*2))
发现这里使用了一个不知道作用的函数 avPlots ，它并没有定义在文件 misc_functions.R 中，可能是R语言的一个函数，不过VSCode并没有将其以红色高亮显示。尝试去掉 avPlots 的作用。先运行第一行代码，看它做了什么：
可知第一行代码从变量fit4.s中取出了PC3.p和wdom_tol这两列数据赋值给变量a。把第一行代码按如下修改，手动取出这两列数据：
a<-data.frame(fit4.s[["model"]][c("PC3.p","wdom_tol")])
names(a)<-c("x","y")
fit4 <- a %>% mutate(y=((scale01(y)-0.5)*2))
使用新的y列数据作为横轴数据再次绘制：
大功告成。不过图像中心附近有几个点的位置并不一致，原因暂时不明。
下面给出“耐湿度(moisture tolerance)”的计算过程：
water01 = Scale01(water.niche.width)
wdom_tol = ranking - water01
moisture_tolerance = ((Scale01(wdom_tol) - 0.5) * 2)
函数 scale01 如下：
scale01<-function(x,na.rm=T,flip=F){
x<-x-min(x,na.rm=T)
x<-(x/max(x,na.rm=T))
if(flip){
x<-(1-x)
}
return(x)
}
即首先减去输入数据的最小值，然后除以此时的最大值，从而把输入数据缩放到 [0,1]。
下面是绘制图C与图A工作的时间表：
2月5日11:03，下载附录A，此后由队友绘制论文A中图C
2月5日15:08，下载论文25及其Github上的数据、代码
2月5日15:12，下载论文34
2月6日09:57，下载并安装R和RStudio，之后安装论文25的代码中用到的R语言包
2月6日11:26，队友提供了绘制出的三张三个温度条件下的图C散点及拟合曲线图
2月6日11:30，论文25的代码运行完毕，并把变量fit3.s导出为几个csv文件
2月6日11:40，导出耐热性数据
2月6日11:44，发现论文25的Github内容新增了两个csv文件 Fungi_moisture_curves.csv 和 Fungi_temperature_curves.csv（论文25的 Fig. 1 c,d 绘图用数据），再次下载（Github显示更新时间为14小时前，大约为2月5日22:00）
2月6日12:00，导出耐湿性数据
2月6日12:32，把三张三个温度条件下的图C散点及拟合曲线图用 StylePix 合成为一张图
2月6日12:40，发现论文25的Github内容更新了 name_list_shortened.csv（新增了两行）文件，再次下载
2月6日13:52，使用队友提供的拟合曲线方程重新绘制，成功画出论文A中图C
2月6日14:08，尝试绘制图A，将三个温度下的分解率分开绘图
2月6日14:33，改用几何平均值绘图
2月6日15:35，导出去掉 avPlots 的作用后的耐湿性数据
2月6日15:50，成功画出论文A中图A
