賣家秀
然而，然而，，我們得到的確實(shí)這樣滴
買家秀
這，這，花紅柳綠的什么玩意?。。?！
大家表示很吃驚，一看文獻(xiàn)原文，是專門挑出了差異基因，難怪這么好看，此時(shí)第一想到的就是測(cè)序公司客服，（微信問）在嗎？好久不見回復(fù)（人家還有大單子，跑你這長(zhǎng)途還不夠油錢）；撥打電話（前提是你有電話號(hào)），沒人接，，，，TMD
于是想自己分析數(shù)據(jù)，然而
linux?
其實(shí)沒有那么難，公司一般都會(huì)提供測(cè)序的count數(shù)據(jù)，自己轉(zhuǎn)化成TPM,FPKM然后挑出顯著差異基因做就O的K啦。
開始展示:
1.獲取基因全長(zhǎng)文件
貼心的我把它放到了網(wǎng)盤，請(qǐng)自取
鏈接：https://pan.baidu.com/s/1FxWPiEwGFobiWka8cNrnbw 提取碼：0cp1 復(fù)制這段內(nèi)容后打開百度網(wǎng)盤手機(jī)App，操作更方便哦--來(lái)自百度網(wǎng)盤超級(jí)會(huì)員V4的分享
2.構(gòu)建帶有對(duì)應(yīng)基因長(zhǎng)度的表達(dá)矩陣
rm(list = ls())#刪除目前工作目錄的變量
library(xlsx)
library(readxl)
ann<- read_excel("Counts.xlsx")#讀取基因文件
input<- read_excel("len.xlsx")#自己在Excel中把網(wǎng)盤里的txt文件基因和長(zhǎng)度共兩列提取出來(lái)
library(dplyr)
merge<-left_join(ann,input,by="Gene")#根據(jù)基因那列進(jìn)行合并
merge <- na.omit(merge)#刪除錯(cuò)誤值行
write.csv(merge,file = "merge.csv",sep = "\t")#讀出文件，直接往下運(yùn)行也許
3,計(jì)算TPM
mycounts<-read.csv("merge.csv")
head(mycounts)
rownames(mycounts)<-mycounts[,1]
mycounts<-mycounts[,-1]
head(mycounts)#最后一列Length是基因長(zhǎng)度
#TPM計(jì)算
kb <- mycounts$Length / 1000
kb
countdata <- mycounts[,1:45]
rpk <- countdata / kb
rpk
tpm <- t(t(rpk)/colSums(rpk) * 1000000)
head(tpm)
write.table(tpm,file="tpm.xls",sep="\t",quote=F)
4.計(jì)算RPKM
fpkm <- t(t(rpk)/colSums(countdata) * 10^6)
head(fpkm)
write.table(fpkm,file="fpkm.xls",sep="\t",quote=F)
Over,是不是super簡(jiǎn)單。
我有一個(gè)愿望，我現(xiàn)在只有人的基因全長(zhǎng)文件，誰(shuí)要是有小鼠的基因全長(zhǎng)文件，一定要聯(lián)系我，我會(huì)把它后續(xù)加上。。。。
好人一生平安，阿門....
