第一天要學會這10題：http://www.bio-info-trainee.com/3750.html
導入文件為Factor時，想轉(zhuǎn)化為character時。這句話不能很少
options(stringsAsFactors = F)
安裝R包就一下幾句命令都試一下就好,包要區(qū)分大小寫
##方法一
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轉(zhuǎn)換尋找到相應的注釋包,下載相應的包
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轉(zhuǎn)換乔询，ID轉(zhuǎn)來轉(zhuǎn)去這三句話就可以（基因ID，探針I(yè)D逢倍，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)
# 根據(jù)名字選擇列
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 數(shù)據(jù)庫挖掘一文就夠
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數(shù)據(jù)分析 網(wǎng)址 http://www.oncolnc.org/
找到TP53基因在TCGA數(shù)據(jù)庫的乳腺癌數(shù)據(jù)集的表達量分組看其是否影響生存
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就這幾句
# 注意查看下載文件的大小扣典，檢查數(shù)據(jù)
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') ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
熱圖如果想要出現(xiàn)臨床信息的注釋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參數(shù)矯正離群值
####apply函數(shù)實戰(zhàn)
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)
繪制相關系數(shù)的熱圖趁怔。（文獻中經(jīng)吃毒耍可將）
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')
對表達矩陣取對數(shù)
expr<- log2(expr+1)
對數(shù)據(jù)進行分組
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的倍數(shù)痕钢，最后2和3的意思是將余數(shù)放到哪里
# 正態(tài)性檢驗
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/胸外/免疫二區(qū)雜志list.xlsx")
list<- list[list$ISSN!= 'NA',]
list<- list[!is.na(list$ISSN) ,] #############
write.csv(list, file = 'magzine_list.csv')
getwd()
image.png
