前言
實現(xiàn)這個功能的軟件也很多置媳，還是煩請大家先自己搜索幾個教程，入門請統(tǒng)一用htseq-count公条，對每個樣本都會輸出一個表達(dá)量文件拇囊。需要用腳本合并所有的樣本為表達(dá)矩陣。參考：生信編程直播第四題：多個同樣的行列式文件合并起來對這個表達(dá)矩陣可以自己簡單在excel或者R里面摸索靶橱，求平均值寥袭，方差」匕裕看看一些生物學(xué)意義特殊的基因表現(xiàn)如何传黄，比如GAPDH,β-ACTIN等等。生信技能樹
htseq-count使用介紹
Usage：htseq-count [options] <sam_file> <gff_file>
主要參數(shù)：
-m 計數(shù)模型選擇队寇，默認(rèn)union膘掰；還有intersection-strict、 intersection-nonempty英上。
-s 是否鏈特異性建庫炭序，默認(rèn)yes
-r 排序方法，推薦按name排序苍日，提高效率
-i 指定read計數(shù)單位惭聂，采用Ensembl GTF時，默認(rèn)按gene_id計數(shù)
-f 輸入文件類型 -a 指定mapping質(zhì)量值過濾read相恃，低的被過濾掉
實驗操作
ls *.Nsort.bam |while read id;
do
# read計數(shù) -s no -r name
nohup samtools view $id |htseq-count -f sam -s no -r name -i gene_name - ../ref/gencode.vM13.annotation.gtf 1>${id%%.*}.geneCounts 2>${id%%.*}.HTseq.log &;
done
# 合并各個樣本的counts值組成表達(dá)矩陣
awk 'BEGIN {OFS="\t"}
NR==FNR{
if(FNR==1){f="Gene\t"FILENAME};
a[FNR]=$1;gene[$1]=$2;next}
{if(FNR==1){f=f"\t"FILENAME};
gene[$1]=gene[$1]"\t"$2}
END{print f;for (i in a) print a[i]"\t"gene[a[i]]}'
*.geneCounts > Akap95_expr.txt
ps: htseq-count目前只能單線程跑辜纲，多個文件通過循環(huán)并行。借助GNU parallel實現(xiàn)htseq-count單文件的并行版拦耐。(存在的問題：個別基因counts數(shù)目比官方的htseq-count多1個count耕腾。原因：對bam分割后，部分reads的兩個mates分割進(jìn)了兩個文件杀糯，htseq-count把這兩個mates都記為了1)
# 借助 GNU parallel扫俺，實現(xiàn) htseq-count 并行版
# -j指定線程；--block文件分塊大小(按1G大小分割)固翰；可根據(jù)具體配置更改
samtools view Akap95_1.Nsort.bam |parallel -j6 --block 1G --pipe -k "htseq-count -f sam -s no -r name -i gene_name - ../ref/gencode.vM13.annotation.gtf >Akap95_1_p{#}.geneCounts" # time ~30min
# 每個job輸出一個文件狼纬，互不干擾
# 用awk合并最終總的counts數(shù)目， !!!個別bam文件分割處的gene骂际，其counts可能有 1 count的差異(相比普通htseq-count)
awk 'BEGIN {OFS="\t"}NR==FNR{a[FNR]=$1;gene[$1]=$2;next}{gene[$1]+=$2} END{for (i in a) print a[i]"\t"gene[a[i]]}' Akap95_1_p*.geneCounts >Akap95_1_p.geneCounts
如下分別是分割后第2個文件最后一行和第3個文件首行：
這兩個reads都比對到了Aggf1的第5號外顯子上
載入IGV看看：(實驗組Akap8(Akap95)表達(dá)更高了疗琉？？歉铝？)
GAPDH盈简、β-actin表達(dá)情況
