215志研代码更新

1+1#1+1代表普通的运算
#install.packages("tidyverse")
#也可以通过右下角的操作栏通过点击来安装R包
#chooseBioCmirror()
#a<-6
#a
#以上是试验赋值的用法，还可以用=来表示，不过==代表一种逻辑判断
#setwd("DAY1")#此代码用来设置工作目录（临时工作目录，退出后回到初始目录），以免在使用过程中将储存文件弄混
#dir.create()可以用来直接创建目录，不过麻烦
#getwd()#查看当前目录
#setwd("DAY1.1")
#ctrl+F可以用来查找代码框中的代码，ctrl+Z是撤销一次的含义
#read.table()#读取文本文档
#write.CSV()#是输出为表格
write.table(counts01A,"counts01A.txt",sep="\t",row.names=T,col.names=NA,quote=F)#是输出文本文档#对文件的读取或输出可以直接点击右上角的读取来导入#通常txt格式会更不容易出错
#t()是用于行列反转的函数
#当一个可以判断为数据框的格式但R语言并没有将其判断为数据框时，可以通过函数as.data.frame()来改变

#class()#用来判断数据框类型，若是数据框内的某一列数据通过输入$符号来确定
#[,]代表对数据框进行提取。其中：表示多少到多少。并且，前是行后是列
#%>%表示传导符号

#install.packages("tidyverse")
#library(tidyverse)
#a<-c("a","b","a","b","c")
#duplicated(a)#用来判断数据里的元素是否重复
#a<-!duplicated(a)#感叹号是指将数据集内判断结果反过来，[]中括号可以把判断中是TRUE的元素提取出来
#a[a]#以上操作中将判断结果赋予了a故导致原本数据的丢失
#a<-c("a","b","a","b","c")
#duplicated(a)
#a<-a[!duplicated(a)]
#a
#这才是正确操作结果

#inner_join()#是将两数据框中数据相同的进行合并的函数，包括（表一，表二，by='合并标准'）
#left_join()#则是代表以函数内左的数据框为标准进行合并


###接下来进行数据的下载实战#####
##建立好文件##
#setwd("TCGA-LUAD")
#setwd('TCGAdata')
library(tidyverse)
install.packages('BiocManager')
library(BiocManager)
BiocManager::install("BioinformaticsFMRP/TCGAbiolinksGUI.data")
BiocManager::install("remotes")
BiocManager::install("ExperimentHub")
BiocManager::install("BioinformaticsFMRP/TCGAbiolinks")
#library(TCGAbiolinks)
cancer_type<-"TCGA-LUAD"
expquery<-GDCquery(project=cancer_type,
                   data.category="Transcriptome Profiling",
                   data.type="Gene Expression Quantification",
                   workflow.type="STAR - Counts")
GDCdownload(expquery,directory="GDCdata")
expquery2<-GDCprepare(expquery,directory="GDCdata",summarizedExperiment = T)
save(expquery2,file="luad.gdc_2023.rad")#对数据进行保存，rad格式是属于R的格式

setwd('TCGA-LUAD')
setwd("TCGAdata")
load("luad.gdc_2023.rad")
load("gene_annotation_2023.rad")
gene_annotation_2023<-data.txt


library(tidyverse)
counts<-expquery2@assays@data@listData[["unstranded"]]
colnames(counts)<-expquery2@colData@rownames
rownames(counts)<-expquery2@rowRanges@ranges@NAMES
counts<-as.data.frame(counts)
counts<-rownames_to_column(counts,var = "ENSEMBL")
counts<-inner_join(counts,gene_annotation_2023,"ENSEMBL")
counts<-counts[!duplicated(counts$symbol),]


rownames(counts)<-NULL#去除行名
counts<-column_to_rownames(counts,var = "symbol")
table(counts$type)#数一下有多少的type
counts<-counts[counts$type == "protein_coding",]
counts<-counts[,-c(1,ncol(counts))]
colnames(counts)<-substring(colnames(counts),1,16)#将列名中1~16字符提取出来 
counts<-counts[,!duplicated(colnames(counts))]



table(substring(colnames(counts),14,16))#通常01A代表肿瘤样本，11A代表正常样本
counts01A<-counts[,substring(colnames(counts),14,16)=="01A"]
counts11A<-counts[,substring(colnames(counts),14,16)=="11A"]#将样本中01A和11A的样本提出来做分析



write.table(counts01A,"counts01A.txt",sep="\t",row.names=T,col.names=NA,quote=F)
write.table(counts11A,"counts11A.txt",sep="\t",row.names=T,col.names=NA,quote=F)
###接下来提取tpms###前面提的counts是用来做差异分析的，而tpms是用来做后面分析的

tpms1<-expquery2@assays@data@listData[["tpm_unstrand"]]
colnames(tpms1)<-expquery2@colData@rownames
rownames(tpms1)<-expquery2@rowRanges@ranges@NAMES
tpms1<-tpms1%>%
  as.data.frame()%>%
  rownames_to_column("ENSEMBL")%>%
  inner_join(gene_annotation_2023,"ENSEMBL")%>%
  .[!duplicated(.$symbol),]
rownames(tpms1)<-NULL
tpms1<-tpms1 %>%column_to_rownames("symbol")
  
table(tpms1$type)#查看分类计数

tpms1<-tpms1[tpms1$type == "protein_coding",]
tpms1<-tpms1[,-c(1,ncol(tpms1))]
colnames(tpms1)<-substring(colnames(tpms1),1,16)#将列名中1~16字符提取出来 
tpms1<-tpms1[,!duplicated(colnames(tpms1))]
table(substring(colnames(tpms1),14,16))

tpms01A<-tpms1[,substring(colnames(tpms1),14,16)=="01A"]
tpms11A<-tpms1[,substring(colnames(tpms1),14,16)=="11A"]


identical(rownames(counts01A),rownames(tpms01A))
identical(rownames(counts11A),rownames(tpms11A))
identical(colnames(counts01A),colnames(tpms01A))
identical(colnames(counts11A),colnames(tpms11A))#验证几个样本是否一致

write.table(tpms01A,"tpms01A.txt",sep="\t",row.names=T,col.names=NA,quote=F)
write.table(tpms11A,"tpms11A.txt",sep="\t",row.names=T,col.names=NA,quote=F)


counts<-cbind(counts01A,counts11A)#cbind是以列colnames来合并，rbind是以行rownames来合并
tpms<-cbind(tpms01A,tpms11A)
write.table(tpms,"tpms.txt",sep="\t",row.names=T,col.names=NA,quote=F)
write.table(counts,"counts.txt",sep="\t",row.names=T,col.names=NA,quote=F)


range(tpms)#查看数据范围
range(tpms01A)
range(tpms11A)
tpms_log2<-log2(tpms+1)
range(tpms_log2)
tpms01A_log2<-log2(tpms01A+1)
range(tpms01A_log2)
tpms11A_log2<-log2(tpms11A+1)
range(tpms11A_log2)
write.table(tpms_log2,"tpms_log2.txt",sep="\t",row.names=T,col.names=NA,quote=F)
write.table(tpms01A_log2,"tpms01A_log2.txt",sep="\t",row.names=T,col.names=NA,quote=F)
write.table(tpms11A_log2,"tpms11A_log2.txt",sep="\t",row.names=T,col.names=NA,quote=F)




###计算肿瘤患者的免疫评分#####
#（