首先列出可能用到的文章：

1. A題題目《2021_MCM_Problem_A》（以下簡(jiǎn)稱題目A）

2. 題目A出題參考《A trait-based understanding of wood decompositionby fungi》（https://www.pnas.org/content/pnas/117/21/11551.full.pdf）（以下簡(jiǎn)稱論文A）

3. 論文A的參考資料25《Consistent trade-offs in fungal trait expression across broad spatial scales》（以下簡(jiǎn)稱論文25）

4. 論文A的參考資料34《Diversity begets diversity in competition for space》（以下簡(jiǎn)稱論文34）
   

然后列出能夠找到的數(shù)據(jù)：

1. 論文A的附錄 (https://www.pnas.org/highwire/filestream/926186/field_highwire_adjunct_files/0/pnas.1909166117.sapp.pdf)（以下簡(jiǎn)稱附錄A）

2. 附錄A的參考資料9 (http://ncf.sobek.ufl.edu/AA00026436)

3. 論文25 "Data availability" 章節(jié)提供的論文所用數(shù)據(jù)和代碼的Github地址 (https://github.com/dsmaynard/fungal_biogeography)

4. 論文34?"Data availability" 章節(jié)提供的論文所用數(shù)據(jù)和代碼的Github地址 (https://github.com/dsmaynard/diversity_begets_diversity)

附錄A提供了S1-S5五張表，其中S3和S4給出了34個(gè)菌種的木材分解率和菌絲延伸率數(shù)據(jù)

使用S3和S4的數(shù)據(jù)就可以畫出題目的圖C了（忽略不確定度）

圖A的繪制就比較麻煩了，需要找到 moisture trade-off 數(shù)據(jù)。

題目A中對(duì)圖A橫軸的說(shuō)明如下：

耐濕性（每種分離物的競(jìng)爭(zhēng)排名及其水分位寬度的差異，均縮放至[0,1]）

moisture tolerance (difference of each isolate’s competitive ranking and their moisture niche width, both scaled to [0,1])

論文A中對(duì)圖A橫軸的說(shuō)明如下：

水分權(quán)衡（A中的x軸）是每個(gè)分離物的競(jìng)爭(zhēng)排名與其水分生態(tài)位寬度之間的差異，兩者均按比例縮放至[0,1]，如Maynard等人(25)所計(jì)算。

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對(duì)左圖(a)的縱軸的說(shuō)明為“耐熱性權(quán)衡” (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提供了數(shù)據(jù)與代碼的下載，其中 Fungal_climate_data.csv 提供了37個(gè)菌種的競(jìng)爭(zhēng)排名 Competitive ranking(ranking)、溫度生態(tài)位寬度 Temperature niche width(temp.niche.width)、水分生態(tài)位寬度 Moisture niche width(water.niche.width) 等數(shù)據(jù)，F(xiàn)ungal_trait_data.csv 除提供上述三種數(shù)據(jù)外，還有生長(zhǎng)率 Growth rate(rate.0.5)?、菌絲密度 Hyphal density(density) 數(shù)據(jù)。需要注意，兩張表格中的競(jìng)爭(zhēng)排名?Competitive ranking(ranking) 數(shù)據(jù)并不一致。

由此推測(cè)，論文作者使用“溫度生態(tài)位寬度”數(shù)據(jù)計(jì)算出“耐熱性”，進(jìn)而繪制出 Fig. 4(a) ，那么改用“水分生態(tài)位寬度”數(shù)據(jù)應(yīng)該就能得到“耐濕性”。為了得知論文作者的計(jì)算過(guò)程，閱讀論文25提供的R語(yǔ)言代碼，可以找到 "#?main?temperature?plot" 注釋的部分和 av_gg 、scale01 等函數(shù)。

通過(guò)網(wǎng)絡(luò)搜索對(duì)代碼中一些R語(yǔ)言關(guān)鍵字的介紹、對(duì)R語(yǔ)言有了簡(jiǎn)單了解后，我安裝了R和RStuido以運(yùn)行代碼。打開(kāi)文件 Regression_and_maps.R ，從頭執(zhí)行至下述語(yǔ)句后，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 的函數(shù)?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數(shù)據(jù)，圖像保存在變量gg2中并作為函數(shù)返回值。利用函數(shù) 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))

運(yùn)行這三行代碼，就得到了保存在變量fit3中的xy數(shù)據(jù)：

保存數(shù)據(jù)到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))

可以發(fā)現(xiàn)，fit3.s與“耐熱性”有關(guān)，fit4.s與“耐濕性”有關(guān)。把之前的三行代碼修改如下并運(yùn)行：

a<-data.frame(avPlots(fit4.s,terms="PC3.p"))

names(a)<-c("x","y")

fit4?<- a %>% mutate(y=((scale01(y)-0.5)*2))

這里的y列數(shù)據(jù)應(yīng)該就是“耐濕度”(moisture tolerance)了。注意到論文25提供的是37個(gè)菌種的數(shù)據(jù)，而附錄A中只有34種，經(jīng)過(guò)對(duì)菌種名稱的比對(duì)，需要去掉編號(hào)16、19、20的菌種數(shù)據(jù)。

觀察圖A，只有一種顏色的點(diǎn)。它的橫軸數(shù)據(jù)應(yīng)該是上面的y列數(shù)據(jù)，但它的縱軸數(shù)據(jù)（對(duì)數(shù)）“分解率”(decomposition rate) 不確定是哪個(gè)溫度條件下的。把三個(gè)溫度下的散點(diǎn)圖都繪制出來(lái)看一下：

觀察最右邊兩個(gè)點(diǎn)，似乎圖A用的是10℃條件下的數(shù)據(jù)，但這兩個(gè)點(diǎn)的縱坐標(biāo)小于2，與圖A不符。查看變量fit4，“耐濕度(moisture tolerance)”為最大值1的菌種編號(hào)36，名稱為 "Tyromyces chioneus HHB11933 B10F"。圖A中該點(diǎn)縱坐標(biāo)約為2.6，e^2.6≈13.46。查看附錄A的S3：

所以圖A的y軸數(shù)據(jù)應(yīng)該是附錄A的S3的最右邊一列“幾何平均值(geometric mean)”。再次繪圖：

還是不對(duì)。再次閱讀代碼：

a<-data.frame(avPlots(fit4.s,terms="PC3.p"))

names(a)<-c("x","y")

fit4?<- a %>% mutate(y=((scale01(y)-0.5)*2))

發(fā)現(xiàn)這里使用了一個(gè)不知道作用的函數(shù)?avPlots ，它并沒(méi)有定義在文件 misc_functions.R 中，可能是R語(yǔ)言的一個(gè)函數(shù)，不過(guò)VSCode并沒(méi)有將其以紅色高亮顯示。嘗試去掉 avPlots 的作用。先運(yùn)行第一行代碼，看它做了什么：

可知第一行代碼從變量fit4.s中取出了PC3.p和wdom_tol這兩列數(shù)據(jù)賦值給變量a。把第一行代碼按如下修改，手動(dòng)取出這兩列數(shù)據(jù)：

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列數(shù)據(jù)作為橫軸數(shù)據(jù)再次繪制：

大功告成。不過(guò)圖像中心附近有幾個(gè)點(diǎn)的位置并不一致，原因暫時(shí)不明。

下面給出“耐濕度(moisture tolerance)”的計(jì)算過(guò)程：

1. water01 = Scale01(water.niche.width)

2. wdom_tol = ranking - water01

3. moisture_tolerance = ((Scale01(wdom_tol) - 0.5) * 2)

函數(shù) 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)

}

即首先減去輸入數(shù)據(jù)的最小值，然后除以此時(shí)的最大值，從而把輸入數(shù)據(jù)縮放到 [0,1]。

下面是繪制圖C與圖A工作的時(shí)間表：

• 2月5日11:03，下載附錄A，此后由隊(duì)友繪制論文A中圖C

• 2月5日15:08，下載論文25及其Github上的數(shù)據(jù)、代碼

• 2月5日15:12，下載論文34

• 2月6日09:57，下載并安裝R和RStudio，之后安裝論文25的代碼中用到的R語(yǔ)言包

• 2月6日11:26，隊(duì)友提供了繪制出的三張三個(gè)溫度條件下的圖C散點(diǎn)及擬合曲線圖
  

• 2月6日11:30，論文25的代碼運(yùn)行完畢，并把變量fit3.s導(dǎo)出為幾個(gè)csv文件

• 2月6日11:40，導(dǎo)出耐熱性數(shù)據(jù)
  

• 2月6日11:44，發(fā)現(xiàn)論文25的Github內(nèi)容新增了兩個(gè)csv文件 Fungi_moisture_curves.csv 和 Fungi_temperature_curves.csv（論文25的 Fig. 1 c,d?繪圖用數(shù)據(jù)），再次下載（Github顯示更新時(shí)間為14小時(shí)前，大約為2月5日22:00）

• 2月6日12:00，導(dǎo)出耐濕性數(shù)據(jù)

• 2月6日12:32，把三張三個(gè)溫度條件下的圖C散點(diǎn)及擬合曲線圖用 StylePix 合成為一張圖

• 2月6日12:40，發(fā)現(xiàn)論文25的Github內(nèi)容更新了 name_list_shortened.csv（新增了兩行）文件，再次下載

• 2月6日13:52，使用隊(duì)友提供的擬合曲線方程重新繪制，成功畫出論文A中圖C
  

• 2月6日14:08，嘗試?yán)L制圖A，將三個(gè)溫度下的分解率分開(kāi)繪圖
  

• 2月6日14:33，改用幾何平均值繪圖

• 2月6日15:35，導(dǎo)出去掉?avPlots 的作用后的耐濕性數(shù)據(jù)

• 2月6日15:50，成功畫出論文A中圖A