Day 1 R语言再复习

第一天要学会这10题：http://www.bio-info-trainee.com/3750.html

导入文件为Factor时，想转化为character时。这句话不能很少

options(stringsAsFactors = F)

安装R包就一下几句命令都试一下就好,包要区分大小写

##方法一
source("http://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
BiocInstaller::biocLite('xCell')  #########这一种方法安装也失败了
########## BioC_mirror: http://bioconductor.org
##方法二
if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
BiocManager::install()
##方法三
install.packages('xCell')
##方法四
install_github("ebecht/MCPcounter",ref="master", subdir="Source")

read.table有很多坑,碰到问题总是在从这里可以解决
csv格式以','分隔

a=read.table('e4-plot.txt',sep = '\t',fill = T,header = T，comment.char = '!', header = T)####真实的表达矩阵
stringsAsFactors = F)####comment.char = '!',    去掉！开头的头文件！

ID转换寻找到相应的注释包,下载相应的包

options()$repos
options()$BioC_mirror 
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
BiocManager::install("请输入自己找到的R包",ask = F,update = F)
options()$repos
options()$BioC_mirror

ID转换，ID转来转去这三句话就可以（基因ID，探针ID，Ensemble ID，Gene Symbol）

library(org.Hs.eg.db)
ls("package:org.Hs.eg.db")
g2s=toTable(org.Hs.egSYMBOL);head(g2s)
g2e=toTable(org.Hs.egENSEMBL);head(g2e)
tmp=merge(a,g2e,by='ensembl_id')
tmp=merge(tmp,g2s,by='gene_id')

文本处理，文字分隔牢记这一句话就够了（使用lapply 或者stringr）

group_list=unlist(lapply(pd$title,function(x){ strsplit(x,' ')[[1]][4]}))
###如果不够再补充几句
####第一句
substring（temprary,1,3）
####第二句stringr 包轻松处理文本
library(stringr)
str_split(patient_ID, ",",simplify = T)[,2] #####stringr直接变成data 

grep筛选行

wax=a[grep('was', a[,1]),]
wax=a[grepl('was', a[,1]),]
was=a[a[,1]=='was',]

dplyr 行列筛选（select 筛选列，filter筛选行，更改名字，%>%管道）

rename(tb, y = year)
filter(iris, Sepal.Length > 7)
# 根据名字选择列
select(flights, year, month, day)

#######"Piping" with %>% makes code more readable, e.g.
iris %>%
group_by(Species) %>%
summarise(avg = mean(Sepal.Width)) %>%
arrange(avg)

data frame 选择行，选择列筛选(match,%in%,merge)三驾马车

######第一架马车
library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
head(ids)
tmp1=merge(ids,a,by='probe_id')
tmp2=ids[match(a$probe_id,ids$probe_id),]

####第二架马车%in%
ng='ACTR3B ANLN BAG1 BCL2 BIRC5 BLVRA CCNB1 CCNE1 CDC20 CDC6'
ng=strsplit(ng,' ')[[1]]
table(ng %in%  rownames(dat))
ng=ng[ng %in%  rownames(dat)]
dat=dat[ng,]
dat=log2(dat)
pheatmap::pheatmap(dat,scale = 'row')

####第三架马车merge
merge
tmp=merge(a,g2e,by='ensembl_id')
tmp=merge(tmp,g2s,by='gene_id')

cbioportal 数据库挖掘一文就够

a=read.table('e4-plot.txt',sep = '\t',fill = T,header = T)
colnames(a)=c('id','subtype','expression','mut')
dat=a
library(ggstatsplot)
ggbetweenstats(data =dat, x = subtype,  y = expression)
library(ggplot2)
ggsave('plot-again-BRCA1-TCGA-BRCA-cbioportal.png')

生存曲线就这两句
单基因生存分析 TCGA数据分析 网址 http://www.oncolnc.org/
找到TP53基因在TCGA数据库的乳腺癌数据集的表达量分组看其是否影响生存

a=read.table('BRCA_7157_50_50.csv',sep = ',',fill = T,header = T)
dat=a
library(ggplot2)
library(survival)
library(survminer) 
table(dat$Status)
dat$Status=ifelse(dat$Status=='Dead',1,0)
sfit <- survfit(Surv(Days, Status)~Group, data=dat)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)
ggsave('survival_TP53_in_BRCA_TCGA.png')

下载GEO就这几句

# 注意查看下载文件的大小，检查数据 
f='GSE24673_eSet.Rdata'
library(GEOquery)
# 这个包需要注意两个配置，一般来说自动化的配置是足够的。
if(!file.exists(f)){
  gset <- getGEO('GSE24673', destdir=".",
                 AnnotGPL = F,     ## 注释文件
                 getGPL = F)       ## 平台文件
  save(gset,file=f)   ## 保存到本地
}
load('GSE24673_eSet.Rdata')  ## 载入数据
class(gset)
length(gset)
class(gset[[1]])

热图如果想要出现临床信息的注释annotation 就要引入tmp的data.frame

a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)
group_list=c('rbc','rbc','rbc','rbn','rbn','rbn','rbc','rbc','rbc','normal','normal')
dat[1:4,1:4]
M=cor(dat)
pheatmap::pheatmap(M)
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)

limma差异表达

exprSet=dat
exprSet[1:4,1:4]
# DEG by limma 
suppressMessages(library(limma)) 
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix<-makeContrasts("progres.-stable",levels = design)
contrast.matrix 
##这个矩阵声明，我们要把progres.组跟stable进行差异分析比较
##step1
fit <- lmFit(exprSet,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix) ##这一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2)  ## default no trend !!!
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput) 
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)

选择重复探针中表达量最大的一个

library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
head(ids)
dat=dat[ids$probe_id,]
dat[1:4,1:4] 
ids$median=apply(dat,1,median)
ids=ids[order(ids$symbol,ids$median,decreasing = T),]
ids=ids[!duplicated(ids$symbol),]
dat=dat[ids$probe_id,]
rownames(dat)=ids$symbol
dat[1:4,1:4]  
dim(dat)

热图画差异表达最明显的50个基因还可以用scale参数矫正离群值

####apply函数实战
apply(e,2,mean) ###2是列的意思
apply(head(e),1,mean)#####1是行的意思

boxplot(apply(dat,1,mad))
pheatmap::pheatmap(dat[order(apply(dat,1,mad), decreasing = T)[1:50],])
pheatmap::pheatmap(dat[order(apply(dat,1,mad), decreasing = F)[1:50],])

####也可以写成这样的
boxplot(apply(dat,1,mad))
cg=order(apply(dat,1,mad), decreasing = T)[1:50]
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(dat)
pheatmap::pheatmap(dat[cg,],annotation_col = tmp)

绘制相关系数的热图。（文献中经常可将）

library(corrplot) 
M <- cor( dat ) 
ac=data.frame(g=group_list)
rownames(ac)=colnames(M) 
pheatmap::pheatmap(M,show_colnames =F,show_rownames = F,
         annotation_col=ac,filename = 'heatmap_cor.png')

主成分分析

  library("FactoMineR")
  library("factoextra") 

  dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)
  fviz_pca_ind(dat.pca,
               geom.ind = c("point", "text"), # show points only (nbut not "text")
               col.ind = dat$group_list, # color by groups
               palette = c("#00AFBB", "#E7B800"),
               addEllipses = TRUE, # Concentration ellipses
               legend.title = "Groups"
  )
  ggsave('all_samples_PCA.png')

对表达矩阵取对数

expr<- log2(expr+1)

对数据进行分组

dt <- read.delim(text =
                   '
pcode name gdp
53 云南 1.2815
54 西藏 0.0921
61 陕西 1.7689
62 甘肃 0.6837
63 青海 0.2301
64 宁夏 0.2752
1 hangzhou 0.2
65 新疆 0.9264', sep=' ', header=T)

dt_1 <- dt %>%
  arrange(desc(gdp)) %>%
  mutate(class=c(rep(c(1,2,3), each=2),2,3)) ######先排序后分组，因为与共是8个不是3的倍数，最后2和3的意思是将余数放到哪里
# 正态性检验
shapiro.test(dt_1$gdp)
# 方差分析
gdp_aov <- aov(gdp~class, dt_1)
summary(gdp_aov)

# Kruskal-Wallis检验
kruskal.test(gdp~class, dt_1)

!is.na()的用法

list <- read_excel("C:/Users/chenyuqiao/Desktop/胸外/免疫二区杂志list.xlsx")
list<- list[list$ISSN!= 'NA',]
list<- list[!is.na(list$ISSN) ,]  #############
write.csv(list, file = 'magzine_list.csv')
getwd()

image.png