任務(wù)
1.實現(xiàn)這個功能的軟件也很多滓玖，還是煩請大家先自己搜索幾個教程坪哄，入門請統(tǒng)一用htseq-count，對每個樣本都會輸出一個表達(dá)量文件。
2.需要用腳本合并所有的樣本為表達(dá)矩陣翩肌。參考：生信編程直播第四題：多個同樣的行列式文件合并起來模暗。
3.對這個表達(dá)矩陣可以自己簡單在excel或者R里面摸索，求平均值念祭，方差兑宇。
看看一些生物學(xué)意義特殊的基因表現(xiàn)如何，比如GAPDH,β-ACTIN等等粱坤。
來源于生信技能樹
下面摘自hoptop隶糕，想看更加詳細(xì)的內(nèi)容htseq的使用方法和計算原理點擊這里
在上篇的比對中，我們需要糾結(jié)是否真的需要比對站玄，如果你只需要知道已知基因的表達(dá)情況枚驻，那么可以選擇alignment-free工具（例如salmon, sailfish)，如果你需要找到noval isoforms蜒什，那么就需要alignment-based工具（如HISAT2, STAR）测秸。到了這一篇的基因（轉(zhuǎn)錄本）定量，需要考慮的因素就更加多了灾常，以至于我不知道如何說清才能理清邏輯霎冯。
定量分為三個水平
- 基因水平(gene-level)
- 轉(zhuǎn)錄本水平(transcript-level)
- 外顯子使用水平(exon-usage-level)。
在基因水平上钞瀑，常用的軟件為HTSeq-count沈撞，featureCounts，BEDTools, Qualimap, Rsubread, GenomicRanges等雕什。以常用的HTSeq-count為例缠俺，這些工具要解決的問題就是根據(jù)read和基因位置的overlap判斷這個read到底是誰家的孩子。值得注意的是不同工具對multimapping reads處理方式也是不同的贷岸，例如HTSeq-count就直接當(dāng)它們不存在壹士。而Qualimpa則是一人一份，平均分配偿警。
對每個基因計數(shù)之后得到的count matrix再后續(xù)的分析中躏救，要注意標(biāo)準(zhǔn)化的問題。如果你要比較同一個樣本(within-sample)不同基因之間的表達(dá)情況螟蒸，你就需要考慮到轉(zhuǎn)錄本長度盒使，因為轉(zhuǎn)錄本越長，那么檢測的片段也會更多七嫌，直接比較等于讓小孩和大人進(jìn)行賽跑少办。如果你是比較不同樣本（across sample）同一個基因的表達(dá)情況，雖然不必在意轉(zhuǎn)錄本長度诵原，但是你要考慮到測序深度（sequence depth)英妓，畢竟測序深度越高挽放，檢測到的概率越大。除了這兩個因素外蔓纠，你還需要考慮GC%所導(dǎo)致的偏差骂维，以及測序儀器的系統(tǒng)偏差。目前對read count標(biāo)準(zhǔn)化的算法有RPKM（SE）, FPKM（PE）贺纲，TPM, TMM等航闺，不同算法之間的差異與換算方法已經(jīng)有文章進(jìn)行整理和吐槽了。但是猴誊，有一些下游分析的軟件會要求是輸入的count matrix是原始數(shù)據(jù)潦刃，未經(jīng)標(biāo)準(zhǔn)化，比如說DESeq2懈叹，這個時候你需要注意你上一步所用軟件會不會進(jìn)行標(biāo)準(zhǔn)化乖杠。
在轉(zhuǎn)錄本水平上，一般常用工具為Cufflinks和它的繼任者StringTie澄成， eXpress胧洒。這些軟件要處理的難題就時轉(zhuǎn)錄本亞型（isoforms）之間通常是有重疊的，當(dāng)二代測序讀長低于轉(zhuǎn)錄本長度時墨状，如何進(jìn)行區(qū)分卫漫？這些工具大多采用的都是expectation maximization（EM）。好在我們有三代測序肾砂。
上述軟件都是alignment-based列赎，目前許多alignment-free軟件，如kallisto, silfish, salmon镐确，能夠省去比對這一步包吝，直接得到read count，在運行效率上更高源葫。不過最近一篇文獻(xiàn)[1]指出這類方法在估計豐度時存在樣本特異性和讀長偏差诗越。
在外顯子使用水平上，其實和基因水平的統(tǒng)計類似息堂。但是值得注意的是為了更好的計數(shù)嚷狞，我們需要提供無重疊的外顯子區(qū)域的gtf文件[2]。用于分析差異外顯子使用的DEXSeq提供了一個Python腳本（dexseq_prepare_annotation.py）執(zhí)行這個任務(wù)储矩。小結(jié)
計數(shù)分為三個水平： gene-level, transcript-level, exon-usage-level
標(biāo)準(zhǔn)化方法： FPKM RPKM TMM TPM
htseq的使用方法和計算原理點擊這里
如何判斷一個 reads 屬于某個基因感耙， htseq-count 提供了 union, intersection_strict,intersection_nonempty 3 種模型褂乍，如圖（大多數(shù)情況下作者推薦用 union 模型）持隧，如果這三種模型還是不和你心意，可以通過htseq-count 自定義模型逃片，方法詳見 A tour through HTSeq 屡拨。
The above figure illustrates the effect of these three modes.png
1.數(shù)據(jù)準(zhǔn)備
輸入為sam格式的文件只酥，如果是paired-end（雙端測序），必須對SAM文件按照reads名稱排序（sort by name呀狼，not position）裂允。對read name或 位置 進(jìn)行排序皆可，通過 -r 可以指定您的數(shù)據(jù)是以什么方式排序的： <order> 可以是 name 或 pos 哥艇， 默認(rèn)是name.命令為
samtools sort -n accepted_hits_unique.bam -o accepted_hits_unique_name_sorted.bam # -n 按read name 排序 绝编，如果不指定則按染色體位置(pos)排序
具體來說：先用samtools先對bam文件排序，再轉(zhuǎn)換為sam(這一步可以不要貌踏，不必進(jìn)行轉(zhuǎn)換）十饥。請看Jimmy文章
# 首先將bam文件按reads名稱進(jìn)行排序（前期是按照默認(rèn)的染色體位置進(jìn)行排序的，所以要重新進(jìn)行排序）祖乳，這里我們主要以小鼠的數(shù)據(jù)為例子逗堵，不進(jìn)行人類的測序數(shù)據(jù)。
$ cd /mnt/f/rna_seq/aligned
#人部分
# -n 按read name 排序 眷昆，如果不指定則按染色體位置排序
$ for ((i=56;i<=58;i++));do samtools sort -n SRR35899${i}.bam -o SRR35899${i}_nsorted.bam;done
#小鼠部分
$ for ((i=59;i<=62;i++));do samtools sort -n SRR35899${i}.bam -o SRR35899${i}_nsorted.bam;done
# 將排序后的bam文件再次轉(zhuǎn)換成sam格式的文件
#這一步可以省略不要
$ for ((i=56;i<=62;i++));do samtools view -h SRR35899${i}_nsorted.bam > SRR35899${i}_count.sam;done
2 reads計數(shù)蜒秤，得到表達(dá)矩陣
數(shù)據(jù)準(zhǔn)備已經(jīng)完成，接下來要使用htseq-count對數(shù)據(jù)進(jìn)行reads 計數(shù)亚斋。
Usage:htseq-count [options] <sam_file> <gff_file>
注：這里最好使用ensembl的基因組注釋文件作媚，小鼠注釋文件下載地址。但是也可以用前面下載過的gencode注釋文件帅刊。
安裝htseq
# 安裝htseq
$ conda install htseq
$ htseq-count --v
#usage: htseq-count [options] alignment_file gff_file
看一下htseq-count幫助文件掂骏，了解下用法
$ htseq-count --h
(base) kelly@DESKTOP-MRA1M1F:/mnt/f/rna_seq/data/matrix$ htseq-count --h
usage: htseq-count [options] alignment_file gff_file
This script takes one or more alignment files in SAM/BAM format and a feature
file in GFF format and calculates for each feature the number of reads mapping
to it. See http://htseq.readthedocs.io/en/master/count.html for details.
positional arguments:
samfilenames Path to the SAM/BAM files containing the mapped reads.
If '-' is selected, read from standard input
featuresfilename Path to the file containing the features
optional arguments:
-h, --help show this help message and exit
-f {sam,bam}, --format {sam,bam}
type of <alignment_file> data, either 'sam' or 'bam'
(default: sam)
-r {pos,name}, --order {pos,name}
'pos' or 'name'. Sorting order of <alignment_file>
(default: name). Paired-end sequencing data must be
sorted either by position or by read name, and the
sorting order must be specified. Ignored for single-
end data.
可以看出，htseq支持SAM和BAM文件厚掷，所以就沒必要把name-sorted BAM文件轉(zhuǎn)成SAM文件弟灼，那樣會用去很多時間和空間。這個通過-f bam解決冒黑。另外雙端測序數(shù)據(jù)必須進(jìn)行排序田绑，看-r選項，即支持染色體位置排序（pos）抡爹，又支持reads name排序掩驱，但name排序會更好。前面已經(jīng)按name排序完成
#存放目錄
$ cd mnt/f/rna_seq/data/reference/annotation/mm10
#下載（很快）
$ wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_mouse/release_M10/gencode.vM10.chr_patch_hapl_scaff.annotation.gtf.gz
#解壓
$ gunzip gencode.vM10.chr_patch_hapl_scaff.annotation.gtf.gz && rm -rf gencode.vM10.chr_patch_hapl_scaff.annotation.gtf.gz
#創(chuàng)建矩陣文件夾
$ cd /mnt/f/rna_seq/data && mkdir -p matrix && cd matrix
#使用htseq-count處理一個小鼠文件
$ htseq-count -r name -f bam /mnt/f/rna_seq/aligned/SRR3589959_nsorted.bam /mnt/f/rna_seq/data/reference/annotation/mm10/gencode.vM10.chr_patch_hapl_scaff.annotation.gtf
還可以下面這樣處理
for ((i=59;i<=62;i++));do htseq-count -r name -f bam /mnt/f/rna_seq/aligned/SRR35899${i}_nsorted.bam /mnt/f/rna_seq/data/reference/annotation/mm10/gencode.vM10.chr_patch_hapl_scaff.annotation.gtf > /mnt/f/rna_seq/data/matrix/SRR35899${i}.count;done
或這樣處理
for i in `seq 59 62`
do
htseq-count -r name -f bam /mnt/f/rna_seq/aligned/SRR35899${i}_nsorted.bam /mnt/f/rna_seq/data/reference/annotation/mm10/gencode.vM10.chr_patch_hapl_scaff.annotation.gtf > /mnt/f/rna_seq/data/matrix/SRR35899${i}.count 2> /mnt/f/rna_seq/data/matrix/SRR35899${i}.log
done
一共處理了4個小鼠的數(shù)據(jù)冬竟，大概5h,結(jié)果如下
(base) kelly@DESKTOP-MRA1M1F:/mnt/f/rna_seq/data/matrix$ ls -al *.count
-rwxrwxrwx 1 root root 1164655 Aug 1 12:35 SRR3589959.count
-rwxrwxrwx 1 root root 1168522 Aug 1 14:09 SRR3589960.count
-rwxrwxrwx 1 root root 1165311 Aug 1 15:17 SRR3589961.count
-rwxrwxrwx 1 root root 1166505 Aug 1 16:39 SRR3589962.count
- 對結(jié)果進(jìn)行統(tǒng)計
$ wc -l *.count
48825 SRR3589959.count
48825 SRR3589960.count
48825 SRR3589961.count
48825 SRR3589962.count
195300 total
- 看下每個文件的格式欧穴，查看前4行,第一列ensembl_gene_id,第二列read_count計數(shù)
$ head -n 4 *.count
==> SRR3589959.count <==
ENSMUSG00000000001.4 1648
ENSMUSG00000000003.15 0
ENSMUSG00000000028.14 835
ENSMUSG00000000031.15 65
==> SRR3589960.count <==
ENSMUSG00000000001.4 2941
ENSMUSG00000000003.15 0
ENSMUSG00000000028.14 1366
ENSMUSG00000000031.15 52
==> SRR3589961.count <==
ENSMUSG00000000001.4 2306
ENSMUSG00000000003.15 0
ENSMUSG00000000028.14 950
ENSMUSG00000000031.15 83
==> SRR3589962.count <==
ENSMUSG00000000001.4 2780
ENSMUSG00000000003.15 0
ENSMUSG00000000028.14 1051
ENSMUSG00000000031.15 53
- 后4行
$ tail -n 4 *.count
==> SRR3589959.count <==
__ambiguous 295498
__too_low_aQual 1107674
__not_aligned 1025857
__alignment_not_unique 11278964
==> SRR3589960.count <==
__ambiguous 498973
__too_low_aQual 1799176
__not_aligned 1675660
__alignment_not_unique 18440618
==> SRR3589961.count <==
__ambiguous 425708
__too_low_aQual 1303141
__not_aligned 1067038
__alignment_not_unique 14080481
==> SRR3589962.count <==
__ambiguous 381237
__too_low_aQual 1542845
__not_aligned 1529309
__alignment_not_unique 14495860
3 合并表達(dá)矩陣并進(jìn)行注釋（R中進(jìn)行）
-
從上面看出需要至少做兩步工作才能更好理解和往下進(jìn)行分析
第一，需要把4個文件合并泵殴；
第二涮帘，需要把ensembl_gene_id轉(zhuǎn)換為gene_symbol;（這一步不進(jìn)行也行，后面還需要）
所以笑诅，上一步得到的4個單獨的矩陣文件调缨，現(xiàn)在要把這4個文件合并為行為基因名疮鲫，列為樣本名，中間為count的矩陣文件弦叶。
先啟動R_studio
(1) 載入數(shù)據(jù)，添加列名
再看下原始數(shù)據(jù),可見59和61和control伤哺，60和62是實驗
> options(stringsAsFactors = FALSE)
> control1<-read.table("SRR3589959.count",sep = "\t",col.names = c("gene_id","control1"))
> head(control1)
gene_id control1
1 ENSMUSG00000000001.4 1648
2 ENSMUSG00000000003.15 0
3 ENSMUSG00000000028.14 835
4 ENSMUSG00000000031.15 65
5 ENSMUSG00000000037.16 70
6 ENSMUSG00000000049.11 0
> control2<-read.table("SRR3589961.count",sep = "\t",col.names = c("gene_id","control2"))
> treat1<-read.table("SRR3589960.count",sep = "\t",col.names = c("gene_id","treat1"))
> treat2<-read.table("SRR3589962.count",sep = "\t",col.names = c("gene_id","treat2"))
(2)數(shù)據(jù)整合
-
merge進(jìn)行整合 -
gencode的注釋文件中的gene_id（如ENSMUSG00000105298.13_3）在EBI是不能搜索到的燕侠，所以用
gsub功能只保留ENSMUSG00000105298這部分
處理之前先看一下,也就是最好5行是我們不需要的，可以刪除
> tail(control1)
gene_id control1
48820 ENSMUSG00000110424.1 26
48821 __no_feature 12642161
48822 __ambiguous 295498
48823 __too_low_aQual 1107674
48824 __not_aligned 1025857
48825 __alignment_not_unique 11278964
當(dāng)然上面這個命令也可以在linux下執(zhí)行
$ tail -n 10 *.count
進(jìn)行整合
> raw_count <- merge(merge(control1, control2, by="gene_id"), merge(treat1, treat2, by="gene_id"))
> head(raw_count)
gene_id control1 control2 treat1 treat2
1 __alignment_not_unique 11278964 14080481 18440618 14495860
2 __ambiguous 295498 425708 498973 381237
3 __no_feature 12642161 15042888 22357626 18675857
4 __not_aligned 1025857 1067038 1675660 1529309
5 __too_low_aQual 1107674 1303141 1799176 1542845
6 ENSMUSG00000000001.4 1648 2306 2941 2780
> tail(raw_count)
gene_id control1 control2 treat1 treat2
48820 ENSMUSG00000110419.1 46 39 68 58
48821 ENSMUSG00000110420.1 0 0 0 0
48822 ENSMUSG00000110421.1 0 0 0 1
48823 ENSMUSG00000110422.1 1 0 1 2
48824 ENSMUSG00000110423.1 0 0 0 0
48825 ENSMUSG00000110424.1 26 20 48 24
刪除前5行
#這里要注意立莉，我讀入之后順序變了贬循，刪除的時候看下刪除的是哪些行
> raw_count_filt <- raw_count[-1:-5,]
gene_id control1 control2 treat1 treat2
6 ENSMUSG00000000001.4 1648 2306 2941 2780
7 ENSMUSG00000000003.15 0 0 0 0
8 ENSMUSG00000000028.14 835 950 1366 1051
9 ENSMUSG00000000031.15 65 83 52 53
10 ENSMUSG00000000037.16 70 53 94 66
11 ENSMUSG00000000049.11 0 3 4 5
因為我們無法在EBI數(shù)據(jù)庫上直接搜索找到ENSMUSG00000024045.5這樣的基因，只能是ENSMUSG00000024045的整數(shù)桃序，沒有小數(shù)點杖虾，所以需要進(jìn)一步替換為整數(shù)的形式。
# 第一步將匹配到的.以及后面的數(shù)字連續(xù)匹配并替換為空媒熊，并賦值給ENSEMBL
>ENSEMBL <- gsub("\\.\\d*", "", raw_count_filt$gene_id)
# 將ENSEMBL重新添加到raw_count_filt1矩陣
>row.names(raw_count_filt) <- ENSEMBL
>head(raw_count_filt)
gene_id control1 control2 treat1 treat2
ENSMUSG00000000001 ENSMUSG00000000001.4 1648 2306 2941 2780
ENSMUSG00000000003 ENSMUSG00000000003.15 0 0 0 0
ENSMUSG00000000028 ENSMUSG00000000028.14 835 950 1366 1051
ENSMUSG00000000031 ENSMUSG00000000031.15 65 83 52 53
ENSMUSG00000000037 ENSMUSG00000000037.16 70 53 94 66
ENSMUSG00000000049 ENSMUSG00000000049.11 0 3 4 5
(3)對基因進(jìn)行注釋-獲取gene_symbol
以下兩種方式可以進(jìn)行
第一：去這里或這里的網(wǎng)頁版奇适，輸入列表即可輸出，不再贅述
第二：用bioMart對ensembl_id轉(zhuǎn)換成gene_symbol
> library('biomaRt')
> library("curl")
Warning message:
package ‘curl’ was built under R version 3.4.4
> mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
> my_ensembl_gene_id<-row.names(raw_count_filt)
> #Sys.setenv(https_proxy="http://proxy:8000")
> options(timeout = 4000000)
> my_ensembl_gene_id<-row.names(raw_count_filt)
> mms_symbols<- getBM(attributes=c('ensembl_gene_id','hgnc_symbol',"chromosome_name", "start_position","end_position", "band"),
+ filters = 'ensembl_gene_id', values = my_ensembl_gene_id, mart = mart)
#于是芦鳍，不幸總是出錯嚷往，我還是用了網(wǎng)頁版工具進(jìn)行轉(zhuǎn)換
> readcount<-read.csv(file="raw_count_filt.csv",header = TRUE)
> head(readcount)
ensembl_id gene_id gene_symbol control1 control2 treat1 treat2
1 ENSMUSG00000000001 ENSMUSG00000000001.4 Gnai3 1648 2306 2941 2780
2 ENSMUSG00000000003 ENSMUSG00000000003.15 Pbsn 0 0 0 0
3 ENSMUSG00000000028 ENSMUSG00000000028.14 Cdc45 835 950 1366 1051
4 ENSMUSG00000000031 ENSMUSG00000000031.15 H19 65 83 52 53
5 ENSMUSG00000000037 ENSMUSG00000000037.16 Scml2 70 53 94 66
6 ENSMUSG00000000049 ENSMUSG00000000049.11 Apoh 0 3 4 5
#查看Akap8(別名AKAP95)表達(dá)情況
> Akap8 <- readcount[readcount$gene_symbol=="Akap8",]
> Akap8
ensembl_id gene_id gene_symbol control1 control2 treat1 treat2
4265 ENSMUSG00000024045 ENSMUSG00000024045.5 Akap8 936 1368 3386 4116
#treat中顯著升高
---------------說明------------------------
以前用bioMart包的getBM命令沒任何問題，這次卻在這里停留了半天了柠衅。一直出下面的錯誤提示皮仁，(你們運行上面代碼很可能不會出錯)
Batch submitting query [===================================>-------------------------------] 53% eta: 4mError in curl::curl_fetch_memory(url, handle = handle) :
Timeout was reached: Connection timed out after 10000 milliseconds
于是不斷google，嘗試解決菲宴，嘗試了以下方法：
- 1 可能是代理問題贷祈，所以采取直接連接和代理連接，未果
- 2 添加
Sys.setenv(https_proxy="http://proxy:8000")喝峦，未果 - 3 添加
options(timeout = 4000000), 未果
至此势誊，我已經(jīng)快崩潰了，不想在這一步一直停留谣蠢，因為這并不影響下一步
又找到一個方法 - 4 把https://www.esemble.org添加到安裝站點粟耻，未果
==================================
后記
經(jīng)過RNA-seq(7): 篩選差異表達(dá)基因（DEGs）部分分析可知，上面代碼沒有問題眉踱，估計是因為我的網(wǎng)實在太不好挤忙，再加上行數(shù)太多（4萬多行）。
==================================
4 輸出count矩陣文件
readcount<-raw_count_filter[ ,-1]
write.csv(readcount, file='readcount.csv')
