作業要求：
實現這個功能的軟件也很多，還是煩請大家先自己搜索幾個教程，入門請統一用htseq-count，對每個樣本都會輸出一個表達量文件。
需要用腳本合並所有的樣本為表達矩陣。參考：生信編程直播第四題：多個同樣的行列式文件合並起來
對這個表達矩陣可以自己簡單在excel或者R里面摸索，求平均值，方差。
看看一些生物學意義特殊的基因表現如何，比如GAPDH,β-ACTIN等等。
【1】安裝計數軟件：htseq-count
1 # conda安裝 2 $ conda install -c bioconda htseq 3 4 # 測試 5 # 能夠在Python中導入HTSeq這個包，說明安裝成功 6 $ python 7 >>> import HTSeq 8 >>>
【2】用htseq對排序后的bam文件進行計數
1 # 用htseq對排序后的bam文件進行計數 2 # htseq-count [options] <alignment_file> <gtf_file> 3 for i in `seq 56 58` 4 do 5 htseq-count -s no -r pos -f bam RNA-Seq/aligned/SRR35899${i}_sorted.bam reference/gencode.v26lift37.annotation.sorted.gtf > RNA-Seq/matrix/SRR35899${i}.count 2> RNA-Seq/matrix/SRR35899${i}.log 6 done 7 8 # （實際情況：沒能對gtf文件進行排序，所以使用未排序的bam文件與原始gtf文件） 9 for i in `seq 56 58` 10 do 11 htseq-count -s no -r pos -f bam hisat2_result/SRR35899${i}.bam reference/genome/gencode/gencode.v26lift37.annotation.gtf > htseq-count_result/SRR35899${i}.count 2>log_SRR35899${i}_htseq 12 done 13 -f bam/sam： 指定輸入文件格式，默認SAM 14 -r name/pos: 你需要利用samtool sort對數據根據read name或者位置進行排序，默認是name 15 -s yes/no/reverse: 數據是否來自於strand-specific assay。DNA是雙鏈的，所以需要判斷到底來自於哪條鏈。如果選擇了no， 那么每一條read都會跟正義鏈和反義鏈進行比較。默認的yes對於雙端測序表示第一個read都在同一個鏈上，第二個read則在另一條鏈上。
程序運行結束后得到矩陣文件：XXXXX.count
【3】合並表達矩陣
1 # R腳本 2 # 1. 導入數據 3 # 首先將四個文件分別賦值：control1，control2，rep1，rep2 4 >control1 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589959.count", sep="\t", col.names = c("gene_id","control1")) 5 >control2 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589961.count", sep="\t", col.names = c("gene_id","control2")) 6 >rep1 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589960.count", sep="\t", col.names = c("gene_id","treat1")) 7 >rep2 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589962.count", sep="\t",col.names = c("gene_id","treat2")) 8 # 2. 數據整合 9 # 將四個矩陣按照gene_id進行合並，並賦值給raw_count 10 >raw_count <- merge(merge(control1, contol2, by="gene_id"), merge(rep1,rep2, by="gene_id")) 11 # 去掉前五行 12 >raw_count_filter <- raw_count[-1:-5, ] 13 # 因為我們無法在EBI數據庫上直接搜索找到ENSMUSG00000024045.5這樣的基因，只能是ENSMUSG00000024045的整數，沒有小數點，所以需要進一步替換為整數的形式。 14 # 第一步將匹配到的.以及后面的數字連續匹配並替換為空，並賦值給ENSEMBL 15 >ENSEMBL <- gsub("\\.\\d*", "", raw_count_filter$gene_id) 16 # 將ENSEMBL重新添加到raw_count_filt1矩陣 17 > row.names(raw_count_filter) <- ENSEMBL 18 # 3. 查看數據整體情況 19 > summary(raw_count_filter) 20 # 4. 將raw_count_filter寫入csv文件（可用Excel文件） 21 > write.csv(raw_count_filter,"d:/cs_file/r_file/file-lyx0623/H_N_rawcount.csv",row.names=FALSE)
理論知識
關於比對：
】想知道已知基因的表達情況：選擇alignment-free工具（例如salmon, sailfish)
】想找noval isoforms：alignment-based工具（如HISAT2, STAR）
】轉錄本定量：需要考慮更多的因素
1、基因表達定量有3個水平：基因水平、轉錄本水平、外顯子使用水平
基因水平：
Htseq-count等工具作用：根據read和基因位置的overlap，判斷該read屬於哪個基因
轉錄本水平：
Cufflinks等工具處理的難題：轉錄本亞型（isoforms）之間通常是有重疊的，當二代測序讀長低於轉錄本長度時，如何進行區分？這些工具大多采用的都是expectation maximization（EM）。好在我們有三代測序。
外顯子使用水平：與基因水平類似。為了更好的計數，需要提供無重疊的外顯子區域的gtf文件，用於分析差異外顯子使用的DEXSeq提供了一個Python腳本（dexseq_prepare_annotation.py）執行這個任務。
2、reads計數后得到的基因定量結果（count matrix），在進行不同維度的比較時需要進行不同的處理（標准化）：
【】比較同一個樣本(within-sample)不同基因之間的表達情況---主要考慮轉錄本長度的影響
（因為轉錄本越長，那么檢測的片段也會更多，直接比較等於讓小孩和大人進行賽跑）
【】比較不同樣本（across-sample）同一個基因的表達情況---主要考慮測序深度的影響
（測序深度越高，檢測到的概率越大）
標准化的算法：RPKM（SE）, FPKM（PE），TPM, TMM，RSEM等等。
一般，標准化之后才能進行差異表達分析（某些軟件要求原始數據，未經標准化）
htseq參數說明
用法：htseq-count [options] alignment_file gff_file
alignment_file：比對結果文件，可以是sam格式或bam格式，可以通過后面介紹的-f參數指定，默認是sam格式。推薦使用按照name排序后的sam文件，這樣的好處是消耗的內存更少，效率更高。關於輸入文件排序詳見參數-r。
gff_file: 包含單位信息的gff/GTF文件（gff文件格式），大多數情況下就是指注釋文件; 由於GTF文件其實就是gff文件格式的變形，在這里同樣可以傳入GTF格式文件。
-f bam/sam： 指定輸入文件格式，默認SAM
-r name/pos: 你需要利用samtool sort對數據根據read name或者位置進行排序，默認是name
-s yes/no/reverse: 數據是否來自於strand-specific assay。DNA是雙鏈的，所以需要判斷到底來自於哪條鏈。如果選擇了no， 那么每一條read都會跟正義鏈和反義鏈進行比較。默認的yes對於雙端測序表示第一個read都在同一個鏈上，第二個read則在另一條鏈上。
-a 最低質量， 剔除低於閾值的read
-m 模式 union（默認）, intersection-strict and intersection-nonempty。一般而言就用默認的，作者也是這樣認為的。
-i id attribute: 在GTF文件的最后一欄里，會有這個基因的多個命名方式（如下）， RNA-Seq數據分析常用的是gene_id， 當然你可以寫一個腳本替換成其他命名方式。