您好,登錄后才能下訂單哦!
本文小編為大家詳細(xì)介紹“R語(yǔ)言怎么實(shí)現(xiàn)數(shù)據(jù)合并和統(tǒng)計(jì)”,內(nèi)容詳細(xì),步驟清晰,細(xì)節(jié)處理妥當(dāng),希望這篇“R語(yǔ)言怎么實(shí)現(xiàn)數(shù)據(jù)合并和統(tǒng)計(jì)”文章能幫助大家解決疑惑,下面跟著小編的思路慢慢深入,一起來(lái)學(xué)習(xí)新知識(shí)吧。
我這做了 一個(gè)對(duì)照和處理各取了四個(gè)時(shí)期的樣品進(jìn)行轉(zhuǎn)錄組分析,每個(gè)樣品3個(gè)生物學(xué)重復(fù),共24個(gè)樣品基因表達(dá)量數(shù)據(jù);
stage1 | stage2 | stage3 | stage4 | ||||
control | DAF2 | DAF5 | DAF11 | DAF16 | |||
case | GDAF2 | GDAF5 | GDAF11 | GDAF16 |
以上是輸入數(shù)據(jù),CK代表DAF,T代表GDAF,樣品24個(gè)有些多,所以換行,后續(xù)分析改了下名字;
1.統(tǒng)計(jì)每個(gè)樣品中表達(dá)的轉(zhuǎn)錄本的個(gè)數(shù),需要對(duì)生物學(xué)重復(fù)樣品基因表達(dá)量進(jìn)行合并(取平均值)
2.繪制柱狀圖,并給出每個(gè)樣品不同范圍基因表達(dá)量柱狀圖
3.繪制樣品之間的PCA分布圖,查看樣品之間的相互關(guān)系,要求不同的處理用相圖的形狀區(qū)分,不同的stage用不同的顏色區(qū)分,相同的stage顏色相同;
4 繪制Venn圖,看看不同的stage 基因的表達(dá)變化等;
1.注意因子的使用,有序因子可以指定柱狀圖樣品的label順序
2.cowplot指定繪圖主題,可直接用于文章發(fā)表主題,SCI主題
3.PCA圖兩種圖例合并,顏色和形狀,方便查看分組和不同的時(shí)期
4.數(shù)據(jù)的合并處理用到apply tapply 非常的方便,避免使用循環(huán)
5 reshape2包的使用,melt把寬型的數(shù)據(jù)轉(zhuǎn)換成了長(zhǎng)型的數(shù)據(jù)方便ggplot2繪圖。
library(reshape2) local({r <- getOption("repos") ;r["CRAN"] <- "http://mirrors.tuna.tsinghua.edu.cn/CRAN/" ;options(repos=r)}) library(ggplot2) library("ggsci") #install.packages('ggedit') library(ggedit) #install.packages("cowplot") library(cowplot) library(Vennerable) library(RColorBrewer) brewer.pal(7,"Set1") setwd("D:/BaiduNetdiskDownload//report/3.gene_expression") getwd() myfpkm<-read.table("All_gene_fpkm.xls",header=TRUE,comment.char="",sep = "\t",check.names=FALSE,row.names=1) head(myfpkm) group=factor(c(rep(c("DAF2", "DAF5", "DAF11", "DAF16"),each=3),rep(c("GDAF2", "GDAF5", "GDAF11", "GDAF16"),each=3)),levels = c("DAF2","DAF5","DAF11","DAF16","GDAF2","GDAF5","GDAF11","GDAF16"),ordered=T) group apply(myfpkm,2,mean) myMeanFun<-function(x){ tapply(as.double(x),group,mean) } meanFpkm=t(apply(myfpkm,1,myMeanFun)) #meanFpkm=meanFpkm[! grepl("newGene",rownames(meanFpkm)),] myCountFun<-function(x){ x=as.double(x) x=x[x>0.5] y=rep("A",times=length(x)) y[x>0.5 & x<=5]<-"0.05<FPKM<=5" y[x>5 & x<100]<-"5<FPKM<100" y[x>=100]<-"FPKM>=100" y table(y) } myCountFPKM=apply(meanFpkm,2,myCountFun) myCountFPKM mycountLongData<-melt(myCountFPKM) pd <- position_dodge(.65) p<-ggplot(mycountLongData)+ geom_bar(mapping = aes(x = factor(Var2,levels = c("DAF2","DAF5","DAF11","DAF16","GDAF2","GDAF5","GDAF11","GDAF16"),ordered=T), y = value, fill=factor(Var1,levels=c("FPKM>=100","5<FPKM<100","0.05<FPKM<=5"),ordered = T)), stat = "identity" ,width=0.5)+ scale_fill_manual(values=c( "#E41A1C", "#377EB8", "#4DAF4A"))+ # theme_bw()+ theme( # panel.grid=element_blank(), # axis.text.x=element_text(colour="black"),axis.text.y=element_text(colour="black"), # panel.border=element_rect(colour = "black"))+ # scale_colour_brewer(palette = 1)+ theme(legend.key = element_blank(),legend.title = element_blank() #axis.text.x = element_text(angle=70, vjust=0.5) )+ xlab("daf")+ylab("count") p pca = prcomp(t(myfpkm[rowSums(myfpkm)>24,]), scale=TRUE) summary(pca) p1=ggplot(as.data.frame(pca$x),aes(x=PC1,y=PC2,colour=group,shape=group))+geom_point(size = 4,alpha=0.7)+ scale_colour_manual(values=c("#E41A1C", "#377EB8", "#4DAF4A" ,"#984EA3","#E41A1C", "#377EB8" ,"#4DAF4A", "#984EA3")) + scale_shape_manual(values=rep(c(17,19),each=4))+ theme(legend.key = element_blank(),legend.title = element_blank()) p1 plot_grid(p, p1, labels = c("A", "B"), align = 'h') ################################################################################### head(meanFpkm) meanFpkm=as.data.frame(meanFpkm) data_list=list(DAF2=geneNames[meanFpkm$DAF2>0.5],DAF5=geneNames[meanFpkm$DAF5>0.5],DAF11=geneNames[meanFpkm$DAF11>0.5],DAF16=geneNames[meanFpkm$DAF16>0.5]) sapply(data_list,length) data<-Venn(data_list) v1=plot(data,doWeight=F,type="ellipses") data_list=list(GDAF2=geneNames[meanFpkm$GDAF2>0.5],GDAF5=geneNames[meanFpkm$GDAF5>0.5],GDAF11=geneNames[meanFpkm$GDAF11>0.5],GDAF16=geneNames[meanFpkm$GDAF16>0.5]) data<-Venn(data_list) v2=plot(data,doWeight=F,type="ellipses")
讀到這里,這篇“R語(yǔ)言怎么實(shí)現(xiàn)數(shù)據(jù)合并和統(tǒng)計(jì)”文章已經(jīng)介紹完畢,想要掌握這篇文章的知識(shí)點(diǎn)還需要大家自己動(dòng)手實(shí)踐使用過(guò)才能領(lǐng)會(huì),如果想了解更多相關(guān)內(nèi)容的文章,歡迎關(guān)注億速云行業(yè)資訊頻道。
免責(zé)聲明:本站發(fā)布的內(nèi)容(圖片、視頻和文字)以原創(chuàng)、轉(zhuǎn)載和分享為主,文章觀點(diǎn)不代表本網(wǎng)站立場(chǎng),如果涉及侵權(quán)請(qǐng)聯(lián)系站長(zhǎng)郵箱:is@yisu.com進(jìn)行舉報(bào),并提供相關(guān)證據(jù),一經(jīng)查實(shí),將立刻刪除涉嫌侵權(quán)內(nèi)容。