題目來(lái)源：http://www.bio-info-trainee.com/3793.html
基于下午的統(tǒng)計(jì)可視化
- 對(duì)前面讀取的
RunInfo Table文件在R里面探索其MBases列头谜，包括 箱線圖(boxplot)和五分位數(shù)(fivenum)蟆盹，還有頻數(shù)圖(hist)，以及密度圖(density) 滥壕。
【相關(guān)知識(shí)鏈接】
箱線圖 - boxplot
五分位數(shù) - fivenum
頻數(shù)圖 - histogram
密度圖 - density
> rm(list = ls())
> options(stringsAsFactors = F)
> sra <- read.table(file = "SraRunTable.txt",header = T,sep = '\t')
#箱線圖
> boxplot(sra$MBases)
箱線圖
#五分位數(shù)
> fivenum(sra$MBases)
[1] 0 8 12 16 74
#頻數(shù)圖
> hist(sra$MBases)
頻數(shù)圖
#密度圖
> density(sra$MBases)
Call:
density.default(x = sra$MBases)
Data: sra$MBases (768 obs.); Bandwidth 'bw' = 1.423
x y
Min. :-4.269 Min. :0.0000000
1st Qu.:16.366 1st Qu.:0.0000353
Median :37.000 Median :0.0003001
Mean :37.000 Mean :0.0121039
3rd Qu.:57.634 3rd Qu.:0.0142453
Max. :78.269 Max. :0.0665647
> plot(density(sra$MBases))
密度圖
#某個(gè)包里畫(huà)密度圖的函數(shù)
> library(lattice)
> densityplot(sra$MBases)
密度圖2
- 把前面讀取的
樣本信息表格的樣本名字根據(jù)下劃線分割看第3列元素的統(tǒng)計(jì)情況变秦。第三列代表該樣本所在的plate 。
> plate <- unlist(strsplit(df_merge$Title,'_'))
> length(plate)
[1] 3072
> tmp <- seq(from = 3,to = length(plate),by = 4)
> plate <- plate[tmp]
> length(plate)
[1] 768
> table(plate)
plate
0048 0049
384 384
根據(jù)plate把關(guān)聯(lián)到的
RunInfo Table信息的MBases列分組檢驗(yàn)是否有統(tǒng)計(jì)學(xué)顯著的差異。
> df_merge$plate <- plate
> t.test(df_merge$MBases~plate)
Welch Two Sample t-test
data: df_merge$MBases by plate
t = 2.3194, df = 728.18, p-value = 0.02065
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.1739531 2.0916719
sample estimates:
mean in group 0048 mean in group 0049
13.78125 12.64844
p<0.05迎膜，有顯著性差異
分組繪制箱線圖(boxplot)，頻數(shù)圖(hist)浆兰，以及密度圖(density) 磕仅。
#箱線圖
>boxplot(df_merge$MBases~df_merge$plate)
箱線圖
#頻數(shù)圖
#按plate拆分df_merge
df_merge_0048 <- df_merge[df_merge$plate=='0048',]
df_merge_0049 <- df_merge[df_merge$plate=='0049',]
hist(df_merge_0048$MBases)
hist(df_merge_0049$MBases)
plate==0048箱線圖
plate==0049箱線圖
#密度圖
plot(density(df_merge_0048$MBases))
plot(density(df_merge_0049$MBases))
plate==0048密度圖
plate==0049密度圖
- 使用ggplot2把上面的圖進(jìn)行重新繪制。
> library(ggplot2)
#箱線圖
ggplot()+
geom_boxplot(
data =df_merge,
mapping = aes(x=plate,y=MBases,fill = plate)
)
#頻數(shù)圖
ggplot()+
stat_count(
data = df_merge,
mapping = aes(x = MBases)
)+
facet_wrap(~plate,scales = "free") #分片
ggplot()+
stat_count(
data = df_merge,
mapping = aes(x = MBases,fill = plate) #填充
)
#密度圖
ggplot()+
geom_density(
data = df_merge,
mapping = aes(x = MBases, fill = plate)
)
ggplot2箱線圖
ggplot2頻數(shù)圖_分片
ggplot2頻數(shù)圖_填充
ggplot2密度圖_分片
ggplot2密度圖_填充
- 使用ggpubr把上面的圖進(jìn)行重新繪制簸呈。
BiocManager::install("ggpubr")
library(ggpubr)
ggboxplot(data = df_merge,x = "plate", y = "MBases",fill = "plate")
gghistogram(data = df_merge,x = "MBases",fill = "plate")
ggdensity(data = df_merge,x = "MBases",fill = "plate")
ggpubr箱線圖
ggpubr頻數(shù)圖
ggpubr密度圖
- 隨機(jī)取384個(gè)MBases信息榕订，跟前面的兩個(gè)plate的信息組合成新的數(shù)據(jù)框，第一列是分組蜕便，第二列是MBases劫恒，總共是384*3行數(shù)據(jù)。
> options(stringsAsFactors = F)
> set.seed(10000)
> id <- sample.int(nrow(df_merge),384,replace = F)
> id
[1] 286 335 301 460 84 265 205 67 575 198 267 272 461 18 422 319 229 226 645 757 558 592 355 520 328 281 442 627 668 641 609 114 152 413 408 3
[37] 39 648 580 135 541 47 295 603 42 14 417 218 707 634 271 638 26 615 162 459 725 154 754 521 390 268 282 386 367 727 637 148 685 765 92 305
[73] 746 658 550 597 312 721 566 686 544 768 58 435 225 248 13 376 684 676 311 85 696 64 481 457 657 759 527 20 533 710 453 10 437 733 595 683
[109] 112 115 735 57 372 583 181 213 697 145 29 705 736 325 208 202 317 687 655 767 652 320 299 425 283 728 485 323 434 734 382 66 186 441 448 196
[145] 432 747 258 339 238 484 375 760 388 297 73 623 83 254 356 654 189 546 32 377 369 37 78 223 157 212 131 147 598 624 451 724 557 633 470 409
[181] 141 547 232 61 694 333 755 287 766 714 589 117 532 76 426 713 499 344 75 748 398 366 555 184 519 664 715 237 309 331 545 48 449 180 296 276
[217] 752 251 74 27 30 581 574 79 610 90 62 682 279 264 310 108 407 54 273 1 554 567 308 36 471 69 672 109 23 496 444 285 81 515 761 304
[253] 509 462 330 630 656 487 34 602 395 132 137 675 418 490 347 242 618 559 616 191 649 552 158 582 190 220 651 53 138 235 172 354 433 342 97 511
[289] 163 401 49 474 763 504 151 531 161 586 421 500 394 628 493 506 65 28 642 660 643 756 368 63 253 525 321 103 24 606 551 40 438 211 577 86
[325] 465 587 19 741 607 72 15 601 133 662 340 343 456 256 56 663 647 362 129 726 447 252 274 599 494 522 187 378 351 718 120 402 743 360 701 611
[361] 244 290 738 45 593 689 518 640 313 125 745 454 165 571 646 17 140 405 680 636 124 25 436 560
> T12 <- data.frame(group = df_merge$plate[id],MBases = df_merge$MBases[id],id = df_merge$Run)
> str(T12)
'data.frame': 768 obs. of 3 variables:
$ group : chr "0048" "0048" "0048" "0049" ...
$ MBases: int 11 15 23 11 17 13 20 21 23 13 ...
$ id : chr "SRR6790720" "SRR6790711" "SRR6790712" "SRR6790713" ...
