轉(zhuǎn)錄組學(xué)習(xí)一（軟件安裝）
轉(zhuǎn)錄組學(xué)習(xí)二（數(shù)據(jù)下載）
轉(zhuǎn)錄組學(xué)習(xí)三（數(shù)據(jù)質(zhì)控）
轉(zhuǎn)錄組學(xué)習(xí)四（參考基因組及gtf注釋探究）
轉(zhuǎn)錄組學(xué)習(xí)五（reads的比對與samtools排序）
轉(zhuǎn)錄組學(xué)習(xí)六（reads計數(shù)與標(biāo)準(zhǔn)化）
轉(zhuǎn)錄組學(xué)習(xí)七（差異基因分析）
轉(zhuǎn)錄組學(xué)習(xí)八（功能富集分析）
任務(wù)
- 學(xué)習(xí)了解各個reads計數(shù)，及標(biāo)準(zhǔn)化的原理，如RPKM/FPKM/TPM的統(tǒng)計學(xué)原理；
- 了解reads計數(shù)的各個軟件，用入門的htseq-count軟件對每個樣本內(nèi)生成關(guān)于表達(dá)量的文件；
- 用腳本合并所有的樣本表達(dá)量文件為表達(dá)矩陣；
- 對表達(dá)矩陣在R軟件里簡單的摸索，例如求平均值，方差等；
- 看一些重要的生物學(xué)意義的特殊基因的表達(dá)情況如何，比如GAPDH,β-ACTIN等。
<font color=orange>reads計數(shù)的原理</font>
參考徐洲更在reads計數(shù)里的解釋：
reads的計數(shù)定量主要可分為三個水平：基因水平、轉(zhuǎn)錄組水平、外顯子水平。
-
<font color=brown>基因水平</font>：常用的軟件包括
HTSeq-count, featureCounts, BEDTools, Qualimap, Rsubread, GenomicRange等。在轉(zhuǎn)錄組學(xué)習(xí)五中的那篇重要的文章對于有參的expression定量所推薦的是FeatureCounts，無參的不組裝直接定量的是Salmon-SMEM。這里對于入門所用的HtSeq-count，解決的就是上一步mapping完得到的比對結(jié)果Sam/Bam文件里的read和有參基因組的overlap區(qū)域來判斷這些reads到底是屬于哪個基因的，然后再進(jìn)一步對這些屬于這個基因的reads進(jìn)行計數(shù)，而出現(xiàn)了overlap區(qū)域就會有不同的情況，對于HTSeq的三種對reads計數(shù)處理的模式如圖：（大多數(shù)情況下，作者推薦union模式）
image
使用完之后，就會對每個樣本輸出一個表達(dá)量的文件，后續(xù)需要對表達(dá)量文件寫腳本合并。 <font color=brown>轉(zhuǎn)錄本水平</font>：使用的工具為Cufflinks和StringTie，eXpress。不同的轉(zhuǎn)錄本之間通常是重疊區(qū)域十分多且相似的，當(dāng)二代讀長小于轉(zhuǎn)錄本長度時，就會有區(qū)分上的難題。不過現(xiàn)在有三代測序，能夠完整的將每個轉(zhuǎn)錄本測得序列。
<font color=brown>外顯子水平</font>：和基因水平定量類似，不過需要提供無重疊的外顯子區(qū)域的gtf文件(?),分析差異外顯子使用的DEXSeq提供了一個Python腳本（dexseq_prepare_annotation.py）執(zhí)行這個任務(wù)
<font color=brown>alignment-free軟件</font>：省去比對，直接得到read count。運(yùn)行效率更高，但會有樣本特異性和讀長偏差（不比對如何知道這些reads是屬于基因組上的哪個gene上，進(jìn)而如何知道定量計數(shù)呢？后續(xù)可以關(guān)注下其算法相關(guān)。）
<font color=orange>標(biāo)準(zhǔn)化的原理RPKM/FPKM/TPM</font>
之前學(xué)習(xí)genek-轉(zhuǎn)錄組原理篇課程時記錄的相關(guān)筆記，也再次回顧一下
-
需要標(biāo)準(zhǔn)化的原因：
- 樣本內(nèi)：相對定量而不是絕對定量：僅表示在35次抽樣中，B基因抽樣到了20次，（而不是其表達(dá)了20次。也不能說其會比geneA表達(dá)量高，因?yàn)榛虻拈L度不同，所以落到基因上的reads數(shù)量也不同）
- 樣本件：不同樣本的測序深度也是不同的，測序深度更深的會得到更多的reads。
- 故，需要標(biāo)準(zhǔn)化。對基因長度和測序深度的不同進(jìn)行標(biāo)準(zhǔn)化。
-
PKM 流程
-
gene長度標(biāo)準(zhǔn)化(真實(shí)轉(zhuǎn)錄序列的長度，不考慮可變剪接情況)
- 測序深度標(biāo)準(zhǔn)化(除以每個樣品總的reads數(shù)目，能夠比對到參考序列上的reads數(shù)目):
-
gene長度標(biāo)準(zhǔn)化(真實(shí)轉(zhuǎn)錄序列的長度，不考慮可變剪接情況)
-
RPKM(Reads Per Kilobase Per Million.)就是上面的對1 million個reads進(jìn)行的操作
- kilobase 基因長度單位
- million 測序深度的單位
- FPKM(Fragments Per Kilobase Per Million).建庫測序是以片段為單位。單端測序兩者等價。而雙端測序應(yīng)該用FPKM是更為通用的表示方式
-
TPM：
- 長度標(biāo)準(zhǔn)化、與FPKM同。
- 測序深度的標(biāo)準(zhǔn)化，(不是如FPKM的除以總reads樹)而是按照長度標(biāo)準(zhǔn)化之后的樣本求和。除以求和的結(jié)果。TPM是相等的
-
gene的表達(dá)量到底代表是什么意思。
- 絕對定量：一個細(xì)胞中，一定mol的RNA種有多少轉(zhuǎn)錄本
- 但是建庫測序時并不知道用于建庫測序時候提取了多少個細(xì)胞或者是總共有多少條轉(zhuǎn)錄本。所以只能進(jìn)行相對定量。
- 相對定量：某一個基因的所有轉(zhuǎn)錄本所占的比例。
-
根據(jù)公式：image
- TPM:總的轉(zhuǎn)錄本的豐度，更符合對相對表達(dá)量的定義。表達(dá)豐度求和。(reads數(shù)目/基因的長度=表達(dá)豐富度)
- 看了許多的文章都在強(qiáng)調(diào)RPKM/FPKM的不合理性跟TPM的相對合理性，但為什么現(xiàn)在大多數(shù)的相對計數(shù)都在用FPKM/RPKM，而不用TPM呢，因?yàn)楹罄m(xù)的軟件不支持。
<font color=orange>HTSeq-count定量</font>
根據(jù)featureCounts or htseq-count?文章里對HTSeq-count和feature-counts做了比較，瀏覽下來得到幾點(diǎn):HTseq-count為Python寫的軟件，為歐洲分子實(shí)驗(yàn)室大神（寫了另外的DESeq2）所寫，并且htseq-count也被更多人所引用。，而對于
featureCounts來說，速度更快，20倍快于htseq；支持SAF的注釋文件，會得到比htseq-count更多的計數(shù)值
- 基本使用：
htseq-count [options] <alignment_file> <gff_file>
- 文獻(xiàn)中測序的數(shù)據(jù)是雙末端PE-reads，htseq的計數(shù)需要進(jìn)行按照reads名稱進(jìn)行排序
### samtools重新排序
for i in `seq 56 58`; do nohup samtools sort -@ 5 -n ../5_samtools/SRR35899${i}.bam -o SRR35899${i}_nsort.bam & done
-
HTSeq-count的基本參數(shù)：
- -f bam/sam:指定文件格式，默認(rèn)為sam，選擇為bam。
- -r name/pos:利用samtool sort對數(shù)據(jù)根據(jù)read name或者位置進(jìn)行排序，默認(rèn)是name。
- -s yes/no/reverse：數(shù)據(jù)是否來自于strand-specific assay。DNA是雙鏈的，所以需要判斷到底來自于哪條鏈。如果選擇了no， 那么每一條read都會跟正義鏈和反義鏈進(jìn)行比較。默認(rèn)的yes對于雙端測序表示第一個read都在同一個鏈上，第二個read則在另一條鏈上。
- -a [num]: 最低質(zhì)量， 剔除低于閾值的read
- -m 模式：union,intersection-strict,intersection-nonempty對應(yīng)著上圖的三種模式。
- **-i **：GFF attribute to be used as feature ID (default,suitable for Ensembl GTF files: gene_id)
- HTSEQ-count腳本：
for i in `seq 56 58`
do
htseq-count -s no -r name -f bam SRR35899${i}_nsort.bam ../../Database/human/Homo_sapiens.GRCh38.87.chr.gtf > SRR35899${i}_matrix.count 2> SRR35899${i}.log
done
-
奇怪的錯誤：
在跑的時候總是存在libbz2.so.1.0不存在，奇奇怪怪的linux文庫問題，在之前安裝samtools=1.6版本時候也遇到這樣情況。此問題有待解決。特么花了兩個小時的時間就在處理這個問題，最后發(fā)現(xiàn)可能是conda軟件在安裝pysam跟htseq時候會有庫的問題，換用imagepip install HTSeq即可解決問題。 -
運(yùn)行結(jié)果：
image
<font color=orange>featureCounts定量</font>
- 基本使用Usage：
featureCounts [options] -a <annotation_file> -o <output_file> input_file1 [input_file2] ...
### 參考其他命令行
featureCounts -T 4 -t CDS -s -g Name -C -a gtf -o readscount_cdf_out 1.sam 3.sam 5.sam ...
$featureCounts -T 5 -p -t exon -g gene_id -a $mm10_gtf -o counts.txt *.bam
<font color=orange>腳本合并各個表達(dá)量的文件生成為表達(dá)矩陣</font>
- Perl里利用哈希數(shù)組來存入數(shù)據(jù)，其中基本知識點(diǎn)邏輯：
- 匹配.count文件的文件名，作為后續(xù)合并的header，
- 將每個基因名做key，value為一維數(shù)組，n->[num];
- 每個hash依次輸出,如果不存在則復(fù)制為"0"
- foreach
#ARGV-1)
-
line[0]}[
line[1];
- print
hash{$id}}),"\n"
#!/usr/bin/perl -w
$header= "Gene";
foreach $i ( 0..@ARGV-1){
#### $ARGV[$i]=~/(.*?)_(.*?).txt/;
$ARGV[$i]=~/(.*?)_matrix\.count/;
$header.="\t$1";
open IN,"$ARGV[$i]" or die "$!";
while(<IN>){
chomp;
@line=split;
$hash{$line[0]}[$i]=$line[1];
}
}
print"$header\n";
foreach $gene_id(sort keys %hash){
foreach $i(0..@ARGV-1){
if(!$hash{$gene_id}[$i]){
$hash{$gene_id}[$i]=0;
}
}
print $gene_id,"\t",join("\t",@{$hash{$gene_id}}),"\n";
}
-
結(jié)果：
image
<font color=orange>重新對小鼠的4個數(shù)據(jù)進(jìn)行比對-計數(shù)</font>
### 下載小鼠的hisat2 參考基因組index數(shù)據(jù)
wget -b ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/grcm38.tar.gz
### 下載小鼠的gtf注釋數(shù)據(jù)：
wget -b wget ftp://ftp.ensembl.org/pub/release-90/gtf/mus_musculus/Mus_musculus.GRCm38.90.gtf.gz
### 對小鼠的數(shù)據(jù)進(jìn)行比對，直接輸出為bam文件，跳過輸出的sam文件。
hisat2 -p 5 -x /sas/supercloud-kong/wangtianpeng/Database/mouse/grcm38/genome -1 /sas/supercloud-kong/wangtianpeng/data/biotrainee/1_raw_data/SRR3589959_1.fastq.gz -2 /sas/supercloud-kong/wangtianpeng/data/biotrainee/1_raw_data/SRR3589959_2.fastq.gz 2>/sas/supercloud-kong/wangtianpeng/data/biotrainee/5_samtools/SRR3589959.log | samtools view -Sb - > SRR3589959.bam
### HISAT2比對的循環(huán)腳本
for i in `seq 60 62`
do
nohup hisat2 -x /sas/supercloud-kong/wangtianpeng/Database/mouse/grcm38/genome -1 /sas/supercloud-kong/wangtianpeng/data/biotrainee/1_raw_data/SRR35899${i}_1.fastq.gz -2 /sas/supercloud-kong/wangtianpeng/data/biotrainee/1_raw_data/SRR35899${i}_2.fastq.gz 2>/sas/supercloud-kong/wangtianpeng/data/biotrainee/5_samtools/SRR35899${i}.log | samtools view -Sb - > /sas/supercloud-kong/wangtianpeng/data/biotrainee/5_samtools/SRR35899${i}.bam &
done
### HTSeq-count計數(shù)
for i in `seq 59 62`
do
nohup htseq-count -s no -r name -f bam -i gene_name ../5_samtools/SRR35899${i}.bam /sas/supercloud-kong/wangtianpeng/Database/mouse/Mus_musculus.GRCm38.90.gtf >SRR35899${i}_matrix.count 2> SRR35899${i}.log &
done
11/22晚11點(diǎn)。Raw_count矩陣導(dǎo)入到R里的探索放在之后的差異基因分析里面吧?；仡欀R點(diǎn)，主要是將HTSeq定量的原理，標(biāo)準(zhǔn)化FPKM/RPKM/TPM的原理，HTSeq-count軟件的使用，還有重新將59~62的小鼠的數(shù)據(jù)跑了一遍?；瞬簧贂r間，好在也學(xué)到了不少。慢慢向前走吧。
參考文檔
- 【使用HTSeq進(jìn)行有參轉(zhuǎn)錄組的表達(dá)量計算】http://www.chenlianfu.com/?p=2438
- 【HOPTOP-reads計數(shù)】https://mp.weixin.qq.com/s?__biz=MzI1MjU5MjMzNA==&mid=2247484502&idx=1&sn=ae6e08bc8ac9374799436f3d301a09ec&chksm=e9e02df7de97a4e14077ca04ac6727ac06012e7d3c3b1ce189bafd679e31b23507078bc0e548&scene=21#wechat_redirect
- 【轉(zhuǎn)錄組入門（6）：reads計數(shù)】http://www.itdecent.cn/p/24cf44b610a7
- 【genek-轉(zhuǎn)錄組原理篇】http://www.genek.tv/
- 【featureCounts or htseq-count?】http://bioinformatics.cvr.ac.uk/blog/featurecounts-or-htseq-count/
- 【轉(zhuǎn)錄組HTseq對基因表達(dá)量進(jìn)行計數(shù)】http://www.bio-info-trainee.com/244.html
