R在數(shù)量生態(tài)學分析中的應用(一)

Numerical Ecology with R 這本書2018年出了第二版黑毅,第一版有中文版的書籍。自己的課題涉及微生物生態(tài)學灶体,但時常對生態(tài)學中的一些專業(yè)名詞概念和分析方法不清楚不理解阅签,嚴重阻礙了對微生物生態(tài)學問題的解決和分析。R語言是生物統(tǒng)計分析和可視化方面的神器蝎抽,感覺這本書在R語言和生態(tài)學的學習中可以起到一舉兩得的作用政钟。

1.數(shù)據(jù)的來源
本書的作者提供了兩個數(shù)據(jù)集以及對應的代碼,對Doubs River流域30個采樣點的魚類以及其他環(huán)境參數(shù)進行收集樟结,形成了Doubs.RData數(shù)據(jù)集养交,對一個區(qū)域環(huán)境70個采樣點中的螨蟲等微小生物以及其他環(huán)境參數(shù)進行收集,形成了mite.RData數(shù)據(jù)集狭吼。

2.數(shù)據(jù)探索
現(xiàn)如今层坠,大多數(shù)生態(tài)學研究需要對數(shù)據(jù)進行假設檢驗和模型化。因此對多維數(shù)據(jù)統(tǒng)計分析的早期是需要使用簡單的統(tǒng)計分析和可視化工具進行數(shù)據(jù)探索性分析的刁笙。以期獲得一下信息:

  • 整個數(shù)據(jù)的全貌概況破花。
  • 對數(shù)據(jù)的變量進行轉(zhuǎn)化和重新編輯
  • 確定進一步分析的方向

2.1數(shù)據(jù)導入和查看

# Load required packages
> library(vegan)
> library(RgoogleMaps)
> library(googleVis)
> library(labdsv)
# Source additional functions that will be used later in this
# Chapter. Our scripts assume that files to be read are in
# the working directory.
> source("panelutils.R")
# Load the data. File Doubs.Rdata is assumed to be
# in the working directory
> load("Doubs.RData")
# The file Doubs.RData contains the following objects:
# spe: species (community) data frame (fish abundances)
# env: environmental data frame
# spa: spatial data frame – cartesian coordinates
# fishtraits: functional traits of fish species
# latlong: spatial data frame – latitude and longitude

## Exploration of a data frame using basic R functions
> spe # Display the whole data frame in the
> spe
   Cogo Satr Phph Babl Thth Teso Chna Pato Lele Sqce Baba Albi Gogo Eslu Pefl Rham Legi Scer Cyca Titi Abbr Icme Gyce
1     0    3    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
2     0    5    4    3    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
3     0    5    5    5    0    0    0    0    0    0    0    0    0    1    0    0    0    0    0    0    0    0    0
4     0    4    5    5    0    0    0    0    0    1    0    0    1    2    2    0    0    0    0    1    0    0    0
5     0    2    3    2    0    0    0    0    5    2    0    0    2    4    4    0    0    2    0    3    0    0    0
6     0    3    4    5    0    0    0    0    1    2    0    0    1    1    1    0    0    0    0    2    0    0    0
7     0    5    4    5    0    0    0    0    1    1    0    0    0    0    0    0    0    0    0    0    0    0    0
8     0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
9     0    0    1    3    0    0    0    0    0    5    0    0    0    0    0    0    0    0    0    1    0    0    0
10    0    1    4    4    0    0    0    0    2    2    0    0    1    0    0    0    0    0    0    0    0    0    0
......
# console
# Not recommended for large datasets!
> spe[1:5, 1:10] # Display only 5 lines and 10 columns
  Cogo Satr Phph Babl Thth Teso Chna Pato Lele Sqce
1    0    3    0    0    0    0    0    0    0    0
2    0    5    4    3    0    0    0    0    0    0
3    0    5    5    5    0    0    0    0    0    0
4    0    4    5    5    0    0    0    0    0    1
5    0    2    3    2    0    0    0    0    5    2
> head(spe) # Display only the first 6 lines
 Cogo Satr Phph Babl Thth Teso Chna Pato Lele Sqce Baba Albi Gogo Eslu Pefl Rham Legi Scer Cyca Titi Abbr Icme Gyce Ruru
1    0    3    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
2    0    5    4    3    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
3    0    5    5    5    0    0    0    0    0    0    0    0    0    1    0    0    0    0    0    0    0    0    0    0
4    0    4    5    5    0    0    0    0    0    1    0    0    1    2    2    0    0    0    0    1    0    0    0    0
5    0    2    3    2    0    0    0    0    5    2    0    0    2    4    4    0    0    2    0    3    0    0    0    5
6    0    3    4    5    0    0    0    0    1    2    0    0    1    1    1    0    0    0    0    2    0    0    0    1
  Blbj Alal Anan
1    0    0    0
2    0    0    0
3    0    0    0
4    0    0    0
5    0    0    0
6    0    0    0

> tail(spe) # Display only the last 6 rows
   Cogo Satr Phph Babl Thth Teso Chna Pato Lele Sqce Baba Albi Gogo Eslu Pefl Rham Legi Scer Cyca Titi Abbr Icme Gyce
25    0    0    0    0    0    0    0    0    1    1    0    0    2    1    0    0    0    1    0    0    0    0    1
26    0    0    0    1    0    0    1    0    1    2    2    1    3    2    1    2    2    1    1    3    2    1    4
27    0    0    0    1    0    0    1    1    2    3    4    1    4    4    1    3    3    1    2    5    3    2    5
28    0    0    0    1    0    0    1    1    2    4    3    1    4    3    2    4    4    2    4    4    3    3    5
29    0    1    1    1    1    1    2    2    3    4    5    3    5    5    4    5    5    2    3    3    4    4    5
30    0    0    0    0    0    0    1    2    3    3    3    5    5    4    5    5    3    5    5    5    5    5    5
   Ruru Blbj Alal Anan
25    1    0    3    0
26    4    2    5    2
27    5    4    5    3
28    5    5    5    4
29    5    4    5    4
30    5    5    5    5

> nrow(spe) # Number of rows (sites)
[1] 30

> ncol(spe) # Number of columns (species)
[1] 27

> dim(spe) # Dimensions of the data frame (rows, columns)
[1] 30 27

> colnames(spe) 
 [1] "Cogo" "Satr" "Phph" "Babl" "Thth" "Teso" "Chna" "Pato" "Lele" "Sqce" "Baba" "Albi" "Gogo" "Eslu" "Pefl" "Rham"
[17] "Legi" "Scer" "Cyca" "Titi" "Abbr" "Icme" "Gyce" "Ruru" "Blbj" "Alal" "Anan"

> rownames(spe) # Row labels (objects = sites)
 [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14" "15" "16" "17" "18" "19" "20" "21" "22" "23"
[24] "24" "25" "26" "27" "28" "29" "30"

> summary(spe) # Descriptive statistics for columns
      Cogo           Satr           Phph            Babl            Thth           Teso             Chna    
 Min.   :0.00   Min.   :0.00   Min.   :0.000   Min.   :0.000   Min.   :0.00   Min.   :0.0000   Min.   :0.0  
 1st Qu.:0.00   1st Qu.:0.00   1st Qu.:0.000   1st Qu.:1.000   1st Qu.:0.00   1st Qu.:0.0000   1st Qu.:0.0  
 Median :0.00   Median :1.00   Median :3.000   Median :2.000   Median :0.00   Median :0.0000   Median :0.0  
 Mean   :0.50   Mean   :1.90   Mean   :2.267   Mean   :2.433   Mean   :0.50   Mean   :0.6333   Mean   :0.6  
 3rd Qu.:0.75   3rd Qu.:3.75   3rd Qu.:4.000   3rd Qu.:4.000   3rd Qu.:0.75   3rd Qu.:0.7500   3rd Qu.:1.0  
 Max.   :3.00   Max.   :5.00   Max.   :5.000   Max.   :5.000   Max.   :4.00   Max.   :5.0000   Max.   :3.0  
      Pato             Lele            Sqce            Baba            Albi          Gogo            Eslu      
 Min.   :0.0000   Min.   :0.000   Min.   :0.000   Min.   :0.000   Min.   :0.0   Min.   :0.000   Min.   :0.000  
 1st Qu.:0.0000   1st Qu.:0.000   1st Qu.:1.000   1st Qu.:0.000   1st Qu.:0.0   1st Qu.:0.000   1st Qu.:0.000  
 Median :0.0000   Median :1.000   Median :2.000   Median :0.000   Median :0.0   Median :1.000   Median :1.000  
 Mean   :0.8667   Mean   :1.433   Mean   :1.867   Mean   :1.433   Mean   :0.9   Mean   :1.833   Mean   :1.333  
 3rd Qu.:2.0000   3rd Qu.:2.000   3rd Qu.:3.000   3rd Qu.:3.000   3rd Qu.:1.0   3rd Qu.:3.750   3rd Qu.:2.000  
 Max.   :4.0000   Max.   :5.000   Max.   :5.000   Max.   :5.000   Max.   :5.0   Max.   :5.000   Max.   :5.000  
      Pefl          Rham          Legi             Scer          Cyca             Titi          Abbr       
 Min.   :0.0   Min.   :0.0   Min.   :0.0000   Min.   :0.0   Min.   :0.0000   Min.   :0.0   Min.   :0.0000  
 1st Qu.:0.0   1st Qu.:0.0   1st Qu.:0.0000   1st Qu.:0.0   1st Qu.:0.0000   1st Qu.:0.0   1st Qu.:0.0000  
 Median :0.5   Median :0.0   Median :0.0000   Median :0.0   Median :0.0000   Median :1.0   Median :0.0000  
 Mean   :1.2   Mean   :1.1   Mean   :0.9667   Mean   :0.7   Mean   :0.8333   Mean   :1.5   Mean   :0.8667  
 3rd Qu.:2.0   3rd Qu.:2.0   3rd Qu.:1.7500   3rd Qu.:1.0   3rd Qu.:1.0000   3rd Qu.:3.0   3rd Qu.:1.0000  
 Max.   :5.0   Max.   :5.0   Max.   :5.0000   Max.   :5.0   Max.   :5.0000   Max.   :5.0   Max.   :5.0000  
      Icme          Gyce            Ruru          Blbj            Alal          Anan     
 Min.   :0.0   Min.   :0.000   Min.   :0.0   Min.   :0.000   Min.   :0.0   Min.   :0.00  
 1st Qu.:0.0   1st Qu.:0.000   1st Qu.:0.0   1st Qu.:0.000   1st Qu.:0.0   1st Qu.:0.00  
 Median :0.0   Median :0.000   Median :1.0   Median :0.000   Median :0.0   Median :0.00  
 Mean   :0.6   Mean   :1.267   Mean   :2.1   Mean   :1.033   Mean   :1.9   Mean   :0.90  
 3rd Qu.:0.0   3rd Qu.:2.000   3rd Qu.:5.0   3rd Qu.:1.750   3rd Qu.:5.0   3rd Qu.:1.75  
 Max.   :5.0   Max.   :5.000   Max.   :5.0   Max.   :5.000   Max.   :5.0   Max.   :5.00  

> range(spe) ## Overall distribution of abundances (dominance codes)
# Minimum and maximum of abundance values in the whole data set
[1] 0 5

> apply(spe, 2, range)
     Cogo Satr Phph Babl Thth Teso Chna Pato Lele Sqce Baba Albi Gogo Eslu Pefl Rham Legi Scer Cyca Titi Abbr Icme
[1,]    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
[2,]    3    5    5    5    4    5    3    4    5    5    5    5    5    5    5    5    5    5    5    5    5    5
     Gyce Ruru Blbj Alal Anan
[1,]    0    0    0    0    0
[2,]    5    5    5    5    5

> apply(spe, 2, range)
     Cogo Satr Phph Babl Thth Teso Chna Pato Lele Sqce Baba Albi Gogo Eslu Pefl Rham Legi Scer Cyca Titi Abbr Icme
[1,]    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
[2,]    3    5    5    5    4    5    3    4    5    5    5    5    5    5    5    5    5    5    5    5    5    5
     Gyce Ruru Blbj Alal Anan
[1,]    0    0    0    0    0
[2,]    5    5    5    5    5

> (ab <- table(unlist(spe)))  # Count the cases for each abundance class

  0   1   2   3   4   5 
435 108  87  62  54  64 

> barplot(ab, las = 1, xlab = "Abundance class", ylab = "Frequency",col = gray(5 : 0 / 5)) # Barplot of the distribution, all species confounded

> sum(spe == 0) # Number of absences
[1] 435

> sum(spe == 0) / (nrow(spe) * ncol(spe)) # Proportion of zeros in the community data set
[1] 0.537037

地理位置參數(shù)的可視化

> plot(spa,
+      asp = 1,
+      type = "n",
+      main = "Site Locations",
+      xlab = "x coordinate (km)",
+      ylab = "y coordinate (km)"
+ )  # Geographic coordinates x and y from the spa data frame
> lines(spa, col = "light blue") # Add a blue line connecting the sites along the Doubs River
> text(spa, row.names(spa), cex = 0.8, col = "red") # Add the site labels
> text(68, 20, "Upstream", cex = 1.2, col = "red") # Add text blocks
> text(15, 35, "Downstream", cex = 1.2, col = "red")

使用Google Maps? map展示結(jié)果

> nom <- latlong$Site
> latlong2 <- paste(latlong$LatitudeN, latlong$LongitudeE, sep = ":")
> df <- data.frame(latlong2, nom, stringsAsFactors = FALSE)
> mymap1 <- gvisMap(df,
+                   locationvar = "latlong2",
+                   tipvar = "nom",
+                   options = list(showTip = TRUE)
+ )
> plot(mymap1)

> par(mfrow = c(2,2)) # Divide the plot window into 4 frames, 2 per row
 # Plot four species
> plot(spa,
+      asp = 1,
+      cex.axis = 0.8,
+      col = "brown",
+      cex = spe$Satr,
+      main = "Brown trout",
+      xlab = "x coordinate (km)",
+      ylab = "y coordinate (km)"
+ ) 
> lines(spa, col = "light blue")
> plot(spa,
+      asp = 1,
+      cex.axis = 0.8,
+      col = "brown",
+      cex = spe$Thth,
+      main = "Grayling",
+      xlab = "x coordinate (km)",
+      ylab = "y coordinate (km)"
+ )
> lines(spa, col = "light blue")
> plot(spa,
+      asp = 1,
+      cex.axis = 0.8,
+      col = "brown",
+      cex = spe$Baba,
+      main = "Barbel",
+      xlab = "x coordinate (km)",
+      ylab = "y coordinate (km)"
+ )
> lines(spa, col = "light blue")
> plot(spa,
+      asp = 1,
+      cex.axis = 0.8,
+      col = "brown",
+      cex = spe$Abbr,
+      main = "Common bream",
+      xlab = "x coordinate (km)",
+      ylab = "y coordinate (km)"
+ )
> lines(spa, col = "light blue")

# Compute the number of sites where each species is present
# To sum by columns, the second argument of apply(), MARGIN,
# is set to 2
> spe.pres <- apply(spe > 0, 2, sum)
> spe.pres 
Cogo Satr Phph Babl Thth Teso Chna Pato Lele Sqce Baba Albi Gogo Eslu Pefl Rham Legi Scer Cyca Titi Abbr Icme Gyce 
   8   17   20   24    8    8   12   11   18   25   14   12   20   18   15   11   13   11   12   17    9    7   12 
Ruru Blbj Alal Anan 
  18   10   14   11 
> sort(spe.pres) #排序
Icme Cogo Thth Teso Abbr Blbj Pato Rham Scer Anan Chna Albi Cyca Gyce Legi Baba Alal Pefl Satr Titi Lele Eslu Ruru 
   7    8    8    8    9   10   11   11   11   11   12   12   12   12   13   14   14   15   17   17   18   18   18 
Phph Gogo Babl Sqce 
  20   20   24   25 

> spe.relf <- 100 * spe.pres/nrow(spe) #計算百分比
> spe.relf 
    Cogo     Satr     Phph     Babl     Thth     Teso     Chna     Pato     Lele     Sqce     Baba     Albi     Gogo 
26.66667 56.66667 66.66667 80.00000 26.66667 26.66667 40.00000 36.66667 60.00000 83.33333 46.66667 40.00000 66.66667 
    Eslu     Pefl     Rham     Legi     Scer     Cyca     Titi     Abbr     Icme     Gyce     Ruru     Blbj     Alal 
60.00000 50.00000 36.66667 43.33333 36.66667 40.00000 56.66667 30.00000 23.33333 40.00000 60.00000 33.33333 46.66667 
    Anan 
36.66667 

> par(mfrow = c(1,2))
> hist(spe.pres,
+      main = "Species Occurrences",
+      right = FALSE,
+      las = 1,
+      xlab = "Number of occurrences",
+      ylab = "Number of species",
+      breaks = seq(0, 30, by = 5),
+      col = "bisque"
+ )
> hist(spe.relf,
+      main = "Species Relative Frequencies",
+      right = FALSE,
+      las = 1,
+      xlab = "Frequency of occurrences (%)",
+      ylab = "Number of species",
+      breaks = seq(0, 100, by = 10),
+      col = "bisque"
+ )

> ## Compare sites: species richness
> # Compute the number of species at each site
> # To sum by rows, the second argument of apply(), MARGIN, is
> # set to 1
> sit.pres <- apply(spe > 0, 1, sum)
> sit.pres
 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 
 1  3  4  8 11 10  5  0  5  6  6  6  6 10 11 17 22 23 23 22 23 22  3  8  8 21 22 22 26 21 
> # Sort the results in increasing order
> sort(sit.pres)
 8  1  2 23  3  7  9 10 11 12 13  4 24 25  6 14  5 15 16 26 30 17 20 22 27 28 18 19 21 29 
 0  1  3  3  4  5  5  6  6  6  6  8  8  8 10 10 11 11 17 21 21 22 22 22 22 22 23 23 23 26 
> par(mfrow = c(1, 2))
> # Plot species richness vs. position of the sites along the river
> plot(sit.pres,type = "s",
+      las = 1,
+      col = "gray",
+      main = "Species Richness vs. \n Upstream-Downstream Gradient",
+      xlab = "Site numbers",
+      ylab = "Species richness"
+ )
> text(sit.pres, row.names(spe), cex = .8, col = "red")
## Compare sites: species richness
# Compute the number of species at each site
# To sum by rows, the second argument of apply(), MARGIN, is
# set to 1
sit.pres <- apply(spe > 0, 1, sum)
# Sort the results in increasing order
sort(sit.pres)
par(mfrow = c(1, 2))
# Plot species richness vs. position of the sites along the river
plot(sit.pres,type = "s",
las = 1,
col = "gray",
main = "Species Richness vs. \n Upstream-Downstream Gradient",
xlab = "Site numbers",
ylab = "Species richness"
)
text(sit.pres, row.names(spe), cex = .8, col = "red")
# Use geographic coordinates to plot a bubble map
plot(spa,
asp = 1,
main = "Map of Species Richness",
pch = 21,
col = "white",
bg = "brown",
cex = 5 * sit.pres / max(sit.pres),
xlab = "x coordinate (km)",
ylab = "y coordinate (km)"
)
lines(spa, col = "light blue")
?著作權歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末,一起剝皮案震驚了整個濱河市疲吸,隨后出現(xiàn)的幾起案子座每,更是在濱河造成了極大的恐慌,老刑警劉巖摘悴,帶你破解...
    沈念sama閱讀 222,000評論 6 515
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件峭梳,死亡現(xiàn)場離奇詭異,居然都是意外死亡蹂喻,警方通過查閱死者的電腦和手機葱椭,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 94,745評論 3 399
  • 文/潘曉璐 我一進店門,熙熙樓的掌柜王于貴愁眉苦臉地迎上來口四,“玉大人孵运,你說我怎么就攤上這事÷剩” “怎么了治笨?”我有些...
    開封第一講書人閱讀 168,561評論 0 360
  • 文/不壞的土叔 我叫張陵,是天一觀的道長赤嚼。 經(jīng)常有香客問我旷赖,道長,這世上最難降的妖魔是什么更卒? 我笑而不...
    開封第一講書人閱讀 59,782評論 1 298
  • 正文 為了忘掉前任等孵,我火速辦了婚禮,結(jié)果婚禮上逞壁,老公的妹妹穿的比我還像新娘流济。我一直安慰自己锐锣,他們只是感情好腌闯,可當我...
    茶點故事閱讀 68,798評論 6 397
  • 文/花漫 我一把揭開白布绳瘟。 她就那樣靜靜地躺著,像睡著了一般姿骏。 火紅的嫁衣襯著肌膚如雪糖声。 梳的紋絲不亂的頭發(fā)上,一...
    開封第一講書人閱讀 52,394評論 1 310
  • 那天分瘦,我揣著相機與錄音蘸泻,去河邊找鬼。 笑死嘲玫,一個胖子當著我的面吹牛悦施,可吹牛的內(nèi)容都是我干的。 我是一名探鬼主播去团,決...
    沈念sama閱讀 40,952評論 3 421
  • 文/蒼蘭香墨 我猛地睜開眼抡诞,長吁一口氣:“原來是場噩夢啊……” “哼!你這毒婦竟也來了土陪?” 一聲冷哼從身側(cè)響起昼汗,我...
    開封第一講書人閱讀 39,852評論 0 276
  • 序言:老撾萬榮一對情侶失蹤,失蹤者是張志新(化名)和其女友劉穎鬼雀,沒想到半個月后顷窒,有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體,經(jīng)...
    沈念sama閱讀 46,409評論 1 318
  • 正文 獨居荒郊野嶺守林人離奇死亡源哩,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 38,483評論 3 341
  • 正文 我和宋清朗相戀三年鞋吉,在試婚紗的時候發(fā)現(xiàn)自己被綠了。 大學時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片励烦。...
    茶點故事閱讀 40,615評論 1 352
  • 序言:一個原本活蹦亂跳的男人離奇死亡谓着,死狀恐怖,靈堂內(nèi)的尸體忽然破棺而出崩侠,到底是詐尸還是另有隱情漆魔,我是刑警寧澤,帶...
    沈念sama閱讀 36,303評論 5 350
  • 正文 年R本政府宣布却音,位于F島的核電站改抡,受9級特大地震影響,放射性物質(zhì)發(fā)生泄漏系瓢。R本人自食惡果不足惜阿纤,卻給世界環(huán)境...
    茶點故事閱讀 41,979評論 3 334
  • 文/蒙蒙 一、第九天 我趴在偏房一處隱蔽的房頂上張望夷陋。 院中可真熱鬧欠拾,春花似錦胰锌、人聲如沸。這莊子的主人今日做“春日...
    開封第一講書人閱讀 32,470評論 0 24
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽。三九已至荆忍,卻和暖如春格带,著一層夾襖步出監(jiān)牢的瞬間,已是汗流浹背刹枉。 一陣腳步聲響...
    開封第一講書人閱讀 33,571評論 1 272
  • 我被黑心中介騙來泰國打工叽唱, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留,地道東北人微宝。 一個月前我還...
    沈念sama閱讀 49,041評論 3 377
  • 正文 我出身青樓棺亭,卻偏偏與公主長得像,于是被迫代替她去往敵國和親蟋软。 傳聞我的和親對象是個殘疾皇子镶摘,可洞房花燭夜當晚...
    茶點故事閱讀 45,630評論 2 359

推薦閱讀更多精彩內(nèi)容