劉小澤寫于2018.8.20
利用htseq-count、featureCounts等定量軟件生成的counts計數文件都是單獨的，后續(xù)進行統計需要將它們匯集到一起
實例測試（單行perl結合R語言）
現在有三個count文件，SRR3589956.count、SRR3589957.count、SRR3589958.count
每一個count文件長這樣
第一步 初步合并各個計數文件
先將兩個文件合并，第一列是樣本名，第二列是基因名，第三列是count結果
perl -lne 'if ($ARGV=~/(.*).count/){print "$1\t$_"}' *.count >matrix.count
三個先合并到一個
第二步 R重塑合并的計數矩陣
a=read.csv('matrix.count',header=F,sep="\t")
colnames(a)=c('sample','gene','count')
library(reshape2)
counts=dcast(a,formula=gene~sample)
write.table(counts, file="join.count",sep="\t",quote=FALSE,row.names=FALSE)
可以直接在命令行運行（前提是安裝了R）
R命令
得到的結果
第三步 探索
- 統計各個樣本count值（最值、中位數）
做個統計
-
按照基因來統計counts的平均數
install.packages("dplyr") library(dplyr) gene_mean = group_by(a,gene)%>%summarize(mean=mean(count)) #另外也可以按照sample來統計 #當然，除了mean平均數，還可以求中值median，最值max、min
按基因來統計count
（沒有實際數據）自己生成測試數據
目的：生成a-g.txt 7個文件，每個文件中第一列為基因名，從gene_1到gene_99，第二列是表達量，1000以內隨機整數
#新建測試文件test.sh
perl -le '{print "gene_$_\t".int(rand(1000)) foreach 1..99}'
# -e后面緊跟著引號里面的字符串是要執(zhí)行的命令；
# -l輸出換行
#想要輸出gene_1這樣類型的，gene_后面的數字用一個變量$_代替，這個變量相當于占一個位置，它的賦值是在后面foreach 1..99，表示 $_從1到99
#加一個\t表示一個tab鍵，空4格
#rand生成隨機數，int限制整數
#注意到"gene_$_\t"與后面生成整數函數之間有一個.【不加.就會報錯，它的作用應該是分離字符串和函數】
#perl腳本調用test.sh
perl -e 'system("bash test.sh > $_.txt") foreach a..g'
#結果就生成了想要的測試文件
開始統計：
#先合并
perl -lne 'if ($ARGV=~/(.*).txt/){print $1\t$_}' *.txt > matrix.count
#進入R studio
#然后重復上面??第二步
歡迎關注我們的公眾號～_～
我們是兩個農轉生信的小碩，打造生信星球，想讓它成為一個不拽術語、通俗易懂的生信知識平臺。需要幫助或提出意見請后臺留言或發(fā)送郵件到Bioplanet520@outlook.com
Welcome to our bioinfoplanet!
