作業(yè)要求
1.實現(xiàn)這個功能的軟件也很多梗逮，還是煩請大家先自己搜索幾個教程秸谢，入門請統(tǒng)一用htseq-count纱新，對每個樣本都會輸出一個表達量文件由捎。
2.需要用腳本合并所有的樣本為表達矩陣商叹。參考：生信編程直播第四題：多個同樣的行列式文件合并起來燕刻。
3.對這個表達矩陣可以自己簡單在excel或者R里面摸索，求平均值剖笙，方差卵洗。
看看一些生物學意義特殊的基因表現(xiàn)如何，比如GAPDH,β-ACTIN等等弥咪。
來源于生信技能樹：http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost
實驗過程
1.數(shù)據(jù)準備
文獻中測序的是PE过蹂，因此我們在對雙端數(shù)據(jù)進行處理時，必須要按reads名進行排序聚至。
# 首先將bam文件按reads名稱進行排序（前期是按照默認的染色體位置進行排序的酷勺，所以要重新進行排序），這里我們主要以小鼠的數(shù)據(jù)為例子扳躬，不進行人類的測序數(shù)據(jù)脆诉。
$ cd ~/disk2/data/rna-seq/aligned
$ 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=59;i<=62;i++));do samtools view -h SRR35899${i}_nsorted.bam > SRR35899${i}_count.sam;done
2. reads計數(shù)
數(shù)據(jù)準備已經(jīng)完成，接下來要使用htseq-count對數(shù)據(jù)進行reads 計數(shù)贷币。
Usage:htseq-count [options] <sam_file> <gff_file>
這里最好使用ensembl的基因組注釋文件击胜，小鼠的文件還是需要再次的下載。
#小鼠基因組注釋文件的下載,我下載的是m10版本的役纹，與基因組匹配
$ mkdir -p ~/disk2/data/reference/genome/mm10
$ cd ~/disk2/data/reference/genome/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
# 之前的環(huán)境已經(jīng)全部搭建完成偶摔，直接就可以使用htseq-count
$ mkdir -p ~/disk2/data/rna-seq/matrix && cd ~/disk2/data/rna-seq/matrix
$ for ((i=59;i<=62;i++));do htseq-count ~/disk2/data/rna-seq/aligned/SRR35899${i}_count.sam ~/disk2/data/reference/genome/mm10/gencode.vM10.chr_patch_hapl_scaff.annotation.gtf;done
這一部分主要參考的是老司機Jimmy大神的博客：http://www.bio-info-trainee.com/244.html
3.合并表達矩陣
使用R將這四個文件進行合并，得到最后的融合表達矩陣字管，R語言代碼：
>options(stringsAsFactors = FALSE)
# 首先將四個文件分別賦值：control1啰挪，control2，rep1嘲叔，rep2
>control1 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589959.count", sep="\t", col.names = c("gene_id","control1"))
>control2 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589961.count", sep="\t", col.names = c("gene_id","control2"))
>rep1 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589960.count", sep="\t", col.names = c("gene_id","akap951"))
>rep2 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589962.count", sep="\t",col.names = c("gene_id","akap952"))
# 將四個矩陣按照gene_id進行合并，并賦值給raw_count
>raw_count <- merge(merge(control1, contol2, by="gene_id"), merge(rep1,rep2, by="gene_id"))
# 需要將合并的raw_count進行過濾處理抽活，里面有5行需要刪除的行硫戈，在我們的小鼠的表達矩陣里面，是1,2,48823,48824,48825這5行下硕。并重新賦值給raw_count_filter
>raw_count_filt <- raw_count[-48823:-48825, ]
>raw_count_filter <- raw_count_filt[-1:-2, ]
# 因為我們無法在EBI數(shù)據(jù)庫上直接搜索找到ENSMUSG00000024045.5這樣的基因丁逝，只能是ENSMUSG00000024045的整數(shù)汁胆，沒有小數(shù)點，所以需要進一步替換為整數(shù)的形式霜幼。
# 第一步將匹配到的.以及后面的數(shù)字連續(xù)匹配并替換為空嫩码，并賦值給ENSEMBL
>ENSEMBL <- gsub("\\.\\d*", "", raw_count_filt1$gene_id)
# 將ENSEMBL重新添加到raw_count_filt1矩陣
>row.names(raw_count_filter) <- ENSEMBL
# 看一些基因的表達情況，在UniProt數(shù)據(jù)庫找到AKAP95的id罪既，并從矩陣中找到訪問铸题，并賦值給AKAP95變量
>AKAP95 <- raw_count_filter[rownames(raw_count_filt1)=="ENSMUSG00000024045",]
# 查看AKAP95
>AKAP95
-
raw_count_filter結(jié)構(gòu)
raw_count_filter數(shù)據(jù)結(jié)構(gòu) - AKAP95基因的表達情況，可以看到表達量是提高
這一部分主要參考徐州更同學的R語言代碼：http://www.reibang.com/p/e9742bbf83b9
我的數(shù)據(jù)是用的小鼠的測序結(jié)果琢感，所以我修改了部分的內(nèi)容丢间，更好的進行。
PS：前段時間一直在忙于實驗驹针，沒有及時去完成這一部分的內(nèi)容烘挫，今天正好是星期天，所以就打算抽出一點時間來完成這一部分的學習筆記柬甥。還有一個原因就是因為自己好像不會了饮六，沒法再進行下去了，說到第還是有些害怕了苛蒲，實在慚愧哈卤橄，繼續(xù)努力加油。
