image.png

前期介绍

癌症类型

全部对比或配对样本；样本数量少时差可从别的数据中进行差异？生存？等后续分析

image.png
样本数量
image.png
image.png

差异分析

counts矩阵必定是整数
image.png
优选count数据，不得已才用以下的数据
image.png

在线数据库

xena

image.png
image.png

数据下载

数据下载

下载方式

image.png

image.png

xena

image.png
image.png
image.png

image.png
image.png

  1. proj = "TCGA-CHOL"
  2. if(F){
  3.   download.file(url = paste0("https://gdc.xenahubs.net/download/",proj, ".htseq_counts.tsv.gz"),destfile = paste0(proj,".htseq_counts.tsv.gz"))
  4.   download.file(url = paste0("https://gdc.xenahubs.net/download/",proj, ".GDC_phenotype.tsv.gz"),destfile = paste0(proj,".GDC_phenotype.tsv.gz"))
  5.   download.file(url = paste0("https://gdc.xenahubs.net/download/",proj, ".survival.tsv"),destfile = paste0(proj,".survival.tsv"))
  6. }

读取和整理数据

查看生存信息与临床信息

生存信息与临床信息在后面的差异分析里用不到，在此仅查看一下，到生存信息部分再整理

  1. clinical = read.delim(paste0(proj,".GDC_phenotype.tsv.gz"),fill = T,header = T,sep = "\t")
  2. surv = read.delim(paste0(proj,".survival.tsv"),header = T)
  3. clinical[1:4,1:4]
  4. #>   submitter_id.samples age_at_initial_pathologic_diagnosis
  5. #> 1     TCGA-ZH-A8Y2-01A                                  59
  6. #> 2     TCGA-ZH-A8Y7-01A                                  59
  7. #> 3     TCGA-W7-A93O-01A                                  NA
  8. #> 4     TCGA-W7-A93O-11A                                  NA
  9. #>   albumin_result_lower_limit albumin_result_specified_value
  10. #> 1                         NA                             NA
  11. #> 2                        3.5                            2.4
  12. #> 3                         NA                             NA
  13. #> 4                         NA                             NA
  14. head(surv)
  15. #>             sample OS    X_PATIENT OS.time
  16. #> 1 TCGA-W5-AA31-01A  0 TCGA-W5-AA31      10
  17. #> 2 TCGA-W5-AA31-11A  0 TCGA-W5-AA31      10
  18. #> 3 TCGA-WD-A7RX-01A  1 TCGA-WD-A7RX      21
  19. #> 4 TCGA-YR-A95A-01A  1 TCGA-YR-A95A      26
  20. #> 5 TCGA-4G-AAZR-11A  1 TCGA-4G-AAZR      28
  21. #> 6 TCGA-4G-AAZR-01A  1 TCGA-4G-AAZR      28

查看矩阵

gdc下载的数据从此处开始衔接

  1. dat = read.table(paste0(proj,".htseq_counts.tsv.gz"),check.names = F,row.names = 1,header = T)
  2. #逆转log
  3. dat = as.matrix(2^dat - 1)
  4. dat[1:4,1:4]
  5. #>                    TCGA-ZD-A8I3-01A TCGA-W5-AA2U-11A TCGA-W5-AA30-01A
  6. #> ENSG00000000003.13             5254             2476             5132
  7. #> ENSG00000000005.5                 1                1                0
  8. #> ENSG00000000419.11             1212              655             1644
  9. #> ENSG00000000457.12              753              346             2652
  10. #>                    TCGA-W5-AA38-01A
  11. #> ENSG00000000003.13             8249
  12. #> ENSG00000000005.5                 1
  13. #> ENSG00000000419.11             1696
  14. #> ENSG00000000457.12              519
  15. # 深坑一个
  16. dat[97,9]
  17. #> [1] 876
  18. as.character(dat[97,9]) #眼见不一定为实吧。
  19. #> [1] "875.999999999999"
  21. # 用apply转换为整数矩阵
  22. exp = apply(dat, 2, as.integer)
  23. exp[1:4,1:4] #行名消失
  24. #>      TCGA-ZD-A8I3-01A TCGA-W5-AA2U-11A TCGA-W5-AA30-01A TCGA-W5-AA38-01A
  25. #> [1,]             5254             2476             5132             8249
  26. #> [2,]                1                1                0                1
  27. #> [3,]             1211              655             1643             1695
  28. #> [4,]              752              345             2652              519
  29. rownames(exp) = rownames(dat)
  30. exp[1:4,1:4]
  31. #>                    TCGA-ZD-A8I3-01A TCGA-W5-AA2U-11A TCGA-W5-AA30-01A
  32. #> ENSG00000000003.13             5254             2476             5132
  33. #> ENSG00000000005.5                 1                1                0
  34. #> ENSG00000000419.11             1211              655             1643
  35. #> ENSG00000000457.12              752              345             2652
  36. #>                    TCGA-W5-AA38-01A
  37. #> ENSG00000000003.13             8249
  38. #> ENSG00000000005.5                 1
  39. #> ENSG00000000419.11             1695
  40. #> ENSG00000000457.12              519
  41. 3.表达矩阵行名ID转换
  42. library(stringr)
  43. head(rownames(exp))
  44. #> [1] "ENSG00000000003.13" "ENSG00000000005.5"  "ENSG00000000419.11"
  45. #> [4] "ENSG00000000457.12" "ENSG00000000460.15" "ENSG00000000938.11"
  46. library(AnnoProbe)
  47. annoGene(rownames(exp),ID_type = "ENSEMBL") # 转换不了
  48. #> Warning in annoGene(rownames(exp), ID_type = "ENSEMBL"): 100% of input IDs are
  49. #> fail to annotate...
  50. #> [1] SYMBOL   biotypes ENSEMBL  chr      start    end     
  51. #> <0 rows> (or 0-length row.names)
  52. rownames(exp) = str_split(rownames(exp),"\\.",simplify = T)[,1];head(rownames(exp))
  53. #> [1] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457"
  54. #> [5] "ENSG00000000460" "ENSG00000000938"
  55. re = annoGene(rownames(exp),ID_type = "ENSEMBL");head(re)
  56. #> Warning in annoGene(rownames(exp), ID_type = "ENSEMBL"): 6.54% of input IDs are
  57. #> fail to annotate...
  58. #>        SYMBOL                           biotypes         ENSEMBL  chr start
  59. #> 1     DDX11L1 transcribed_unprocessed_pseudogene ENSG00000223972 chr1 11869
  60. #> 2      WASH7P             unprocessed_pseudogene ENSG00000227232 chr1 14404
  61. #> 3   MIR6859-1                              miRNA ENSG00000278267 chr1 17369
  62. #> 4 MIR1302-2HG                             lncRNA ENSG00000243485 chr1 29554
  63. #> 6     FAM138A                             lncRNA ENSG00000237613 chr1 34554
  64. #> 7      OR4G4P             unprocessed_pseudogene ENSG00000268020 chr1 52473
  65. #>     end
  66. #> 1 14409
  67. #> 2 29570
  68. #> 3 17436
  69. #> 4 31109
  70. #> 6 36081
  71. #> 7 53312
  72. library(tinyarray)
  73. exp = trans_array(exp,ids = re,from = "ENSEMBL",to = "SYMBOL")
  74. exp[1:4,1:4]
  75. #>             TCGA-ZD-A8I3-01A TCGA-W5-AA2U-11A TCGA-W5-AA30-01A TCGA-W5-AA38-01A
  76. #> DDX11L1                    0                0                0                1
  77. #> WASH7P                    81               10              146               54
  78. #> MIR6859-1                  1                0               11                1
  79. #> MIR1302-2HG                0                0                0                0

image.png
image.png

基因过滤

  1. nrow(exp)#过滤之前基因数
  2. #> [1] 56514
  3. #常用过滤标准1：仅去除在所有样本里表达量都为零的基因
  4. exp1 = exp[rowSums(exp)>0,]
  5. nrow(exp1)
  6. #> [1] 48057
  7. #常用过滤标准2(推荐)：仅保留在一半以上样本里表达的基因
  8. exp = exp[apply(exp, 1, function(x) sum(x > 0) >0.5*ncol(exp)), ]
  9. nrow(exp)
  10. #> [1] 28434

分组信息获取

  1. table(str_sub(colnames(exp),14,15))
  2. #> 
  3. #> 01 11 
  4. #> 36  9
  5. Group = ifelse(as.numeric(str_sub(colnames(exp),14,15)) < 10,'tumor','normal')
  6. Group = factor(Group,levels = c("normal","tumor"))
  7. table(Group)
  8. #> Group
  9. #> normal  tumor 
  10. #>      9     36

保存数据

  1. save(exp,Group,proj,clinical,surv,file = paste0(proj,".Rdata"))

回顾代码

image.png

  1. x = round(rnorm(10),2)
  2. x
  3. sort(x,decreasing = T)
  4. #生成一个与x内容相同顺序不同的向量
  5. y = sample(x,10)
  6. x
  7. match(x,y) #谁的下标
  8. x[match(x,y)] #错
  9. y[match(x,y)]==x
  10. identical(y[match(x,y)],x) # 两个数据是否一致
  12. names(x)=letters[1:10]
  13. class(x)
  15. library(stringr)
  16. x1 = str_split(sentences[2]," ")[[1]]
  17. ifelse(str_detect(x1,"th"),
  18.        "hy",
  19.        "no")
  20. library(tibble)
  21. df1 = as.data.frame(x) %>% 
  22.   rownames_to_column(var = "letters")
  23. library(dplyr)
  24. arrange(df1,desc(x))
  26. df2 = data.frame(letters = letters,
  27.                  y = rnorm(26))
  28. df2
  30. merge(df1,df2,by = "letters")
  32. a = sample(letters[1:5],7,replace = T)
  33. a
  34. unique(a)
  35. a[!duplicated(a)]==unique(a)
  36. identical(a[!duplicated(a)],unique(a))
  38. x[1]
  39. x[x>1]
  40. x["a"]
  42. df1$y = iris$Sepal.Length[1:10]
  43. df1 = mutate(df1,z = iris$Sepal.Width[1:10])
  45. b = read.delim("TCGA-CHOL.GDC_phenotype.tsv.gz")
  47. load("gtf_gene.Rdata")
  48. save(b,file = "b.Rdata")
  49. pdf("a.pdf")
  50. pheatmap::pheatmap(volcano)
  51. dev.off()
  52. library(eoffice)
  53. f = pheatmap::pheatmap(volcano)
  54. topptx(f,"f.pptx")
  55. install.packages()
  56. BiocManager::install()
  57. devtools::install_github()
  58. library()

image.png

image.png
实现仅保留在一半以上样本里表达的基因的过滤标准2(推荐)

  1. exp = exp[apply(exp, 1, function(x) sum(x > 0) > 0.5*ncol(exp)), ] #矩阵按行取子集，表达矩阵筛选基因
  2. nrow(exp)
  3. #对上述代码的理解，如下：
  5. #sum(逻辑值)代表T的个数，T是表达量大于0
  6. #sum(x>0)是某基因在多少个样本里有表达
  7. sum(x1 > 0) > 22.5
  8. sum(x2 > 0) > 22.5
  9. jianan = function(x) sum(x > 0) > 0.5*ncol(exp)
  10. x1 = as.numeric(exp[1,])
  11. x2 = as.numeric(exp[2,])
  12. jianan(x1)
  13. jianan(x2)
  14. #对矩阵的每一行(基因)做了一个函数的操作
  15. k = apply(exp, 1, jianan)
  16. #矩阵按行取子集，表达矩阵筛选基因
  17. k= apply(exp, 1, function(x) sum(x > 0) > 0.5*ncol(exp));table(k)
  18. exp = exp[k,]
  19. dim(exp)
  20. exp = exp[apply(exp, 1, function(x) sum(x > 0) > 0.5*ncol(exp)), ]

image.png
image.png

gdcRNAtools

gdc-client

image.png

image.png
image.png

R markdown

渲染和导出
image.png
image.png
image.png
image.png