2.軟件以及流程以及代碼： 本實(shí)驗(yàn)使用常規(guī)的Tophat2 比對(duì)
（1）樣本制備、建庫(kù)：總RNA ---> polyA富集mRNA ---> 打斷 ---> 隨機(jī)引物反轉(zhuǎn)錄成cDNA ---> 末端修復(fù)、加A寡喝、加接頭
RNA提绕釉颉：方法-使用trizol-based方法或者試劑盒方法聪舒；
富集mRNA掖肋，除去rRNA等RNA:依據(jù)為-在測(cè)mRNA過(guò)程中吐咳，首先要去除rRNA会宪。以人為例肖卧，在抽提的總RNA中，95%的RNA是rRNA掸鹅，2%的RNA是mRNA塞帐，剩下的則是lncRNA、microRNA巍沙、siRNA等葵姥。rRNA整個(gè)人類(lèi)當(dāng)中是非常保守的，在各個(gè)組織器官中也是非常穩(wěn)定的句携，因此這些測(cè)序結(jié)果對(duì)我們的研究是沒(méi)有用處的榔幸。mRNA則是RNA中比較重要的部分。
具體：以人的為例：總RNA跟帶有Poly(T)探針的磁珠結(jié)合-洗脫結(jié)合的mRNA-用Mg離子溶液打碎mRNA-隨機(jī)引物反轉(zhuǎn)第一條鏈的cDNA，之后再合成第二條cDNA,獲得雙鏈cDNA-對(duì)雙鏈cDNA末端修復(fù)削咆，加A加接頭-片段選擇牍疏，PCR擴(kuò)增、純化（如果樣本中存在污染物拨齐，則需要結(jié)合試劑盒進(jìn)一步純化）鳞陨。
###一般在會(huì)測(cè)序前對(duì)總RNA進(jìn)行一次質(zhì)檢，根據(jù)電泳質(zhì)檢結(jié)果中的18S和28S（rRNA）兩個(gè)峰的高度以及峰的尖度來(lái)判斷RNA的質(zhì)量奏黑，峰越高越尖（RIN > 8.0）表示RNA的完整度越好炊邦。當(dāng)然，濃度以及A260/A280比值也是需要的熟史。
?1馁害、真核生物種常規(guī)的去除rRNA的方法是通過(guò)oligo(dT)富集帶有polyA尾的mRNA來(lái)實(shí)現(xiàn)的，2蹂匹、不含有polyA尾的轉(zhuǎn)錄本序列以及存在部分降解的總RNA樣本碘菜，所以這種方法針對(duì)福爾馬林(Formalin-Fixed)樣本和FFPE(Paraffin-Embedded)石蠟包埋樣本是不適用的，否則對(duì)獲得樣本中最全面的轉(zhuǎn)錄本信息會(huì)產(chǎn)生顯著影響限寞，一般采用需結(jié)合RiboZero忍啸、RiboMinus等是結(jié)合來(lái)開(kāi)展去除。針對(duì)FFPE樣本還有結(jié)合雙鏈特異性核酸酶構(gòu)建文庫(kù)來(lái)降低后續(xù)測(cè)序數(shù)據(jù)中的rRNA序列比例的履植。
建庫(kù)：去除rRNA之后獲得的mRNA進(jìn)行構(gòu)建文庫(kù)计雌，先對(duì)mRNA打斷再進(jìn)行反轉(zhuǎn)錄的文庫(kù)構(gòu)建方法。之后反轉(zhuǎn)的cDNA再末端修復(fù)到平末端玫霎，加上ployA和接頭凿滤。
###當(dāng)然，里面涉及的蛋白相關(guān)知識(shí)庶近，例如蛋白變性失活翁脆、鍵的破壞、復(fù)性和鹽析等涉及的蛋白四級(jí)結(jié)構(gòu)鼻种、各種化學(xué)鍵和相互作用的內(nèi)容反番。
（2）測(cè)序：SE測(cè)序-單端測(cè)序；PE測(cè)序-雙端測(cè)序叉钥；一般現(xiàn)在用的是PE測(cè)序罢缸。
##里面有測(cè)序的方法、測(cè)序原理投队、不同品牌的區(qū)別
?.sra-> .fastq? 代碼：fastq-dump ?--gzip ?--split-3 –O ../fastq/ -A ../xx.sra
（3）質(zhì)量控制:當(dāng)然用IGV也可以做一些質(zhì)控枫疆。
Fastqc??自己命名.fastq.gz?? -o保存文件夾/
（4）數(shù)據(jù)與處理：
比對(duì)質(zhì)量過(guò)濾、修剪：trimmomatic? PE輸入文件.1.fastq.gz 輸入文件.2.fastq.gz? paired1.fq.gz unpaired.1.fq.gz? paired2.fq.gz unpaired.2.fq.gz AVGQUAL:20 MINLEN:50? (刪去質(zhì)量小于20蛾洛，且刪除讀段小于50的片段)
除去那些測(cè)不出來(lái)或者說(shuō)未被識(shí)別的序列 標(biāo)記為N，很多的話配對(duì)時(shí)需要?jiǎng)h除這樣的序列：prinseq-lite.pl –fastq read1.fastq
–fastq2 read2.fastq –ns_max_n 2 –out_god nfiltered –out_bad null –no_qual_header
– log –verbose (刪除每條序列上含有2個(gè)N的序列)
去接頭：trimmomatic? PE輸入文件.1.fastq.gz 輸入文件.2.fastq.gz? paired1.fq.gz unpaired.1.fq.gz? paired2.fq.gz unpaired.2.fq.gzILLUMINACLIP:TruSeq2 –PE .fa:2:30:10:1:true?(刪除TruSeq2接頭，允許有2個(gè)不匹配轧膘，回文剪接閾值是30钞螟，簡(jiǎn)單剪接閾值10，回文模式檢測(cè)到的最低接頭長(zhǎng)度是1谎碍，反向讀段被保留-默認(rèn)它被刪除)
重復(fù)：prinseq-lite.pl–fastq read1.fastq –fastq2 read2.fastq –derep_min 101
–out_god nfiltered –out_bad
null –no_qual_header – log –verbose
（5）比對(duì)or 的de novo：
tophat2 –o 保存文件夾/ --transcriptome-index ?轉(zhuǎn)錄本索引bt2/ -p 8 基因組索引bt2/ ?輸入文件
（6）比對(duì)結(jié)果注釋?zhuān)炕疪SeQc
???? 1.比對(duì)統(tǒng)計(jì)項(xiàng)：bam_stat.py –i 比對(duì)文件
???? 2.比對(duì)到基因組各個(gè)原件上的情況：
??????? Read_distribution.py –r 基因組的gtf對(duì)應(yīng)的bed文件–I 比對(duì)后的文件bam
?3.轉(zhuǎn)錄本覆蓋度是不是有偏差
???? geneBody-coverage.py –r 基因組的gtf對(duì)應(yīng)的bed文件–I 比對(duì)后的文件bam
?4.測(cè)序深度上表達(dá)豐度檢測(cè)
? RPKM_saturation.py? –r 基因組的gtf對(duì)應(yīng)的bed文件–I 比對(duì)后的文件bam
5. 剪接點(diǎn)
Junction_annotation.py–r 基因組的gtf對(duì)應(yīng)的bed文件–I 比對(duì)后的文件bam
6.剪接點(diǎn)飽和狀態(tài)檢測(cè)
Junction_saturation.py –r 基因組的gtf對(duì)應(yīng)的bed文件–I 比對(duì)后的文件bam
量化：
1. ? 基因水平上-reads落在那些基因上? HTSeq軟件
比對(duì)結(jié)果按照名字排序：samtools sort –n 比對(duì)文件 排序后的文件
量化：htseq-count –f bam –stranded=no 排過(guò)序的文件bam基因組索引文件gtf > counts.txt?
2. ? 轉(zhuǎn)錄水平量化：cufflinks –G 基因組文件gtf? -b基因組–u –p 8 比對(duì)后為排序的文件bam? -o保存文件夾
3. ? 外顯子水平? DEXSeq軟件 ?
扁平化（即鳞滨，先將基因組注釋文件扁平化，拉開(kāi)距離蟆淀，形成不重疊的外顯子區(qū)域拯啦，進(jìn)而將比對(duì)數(shù)據(jù)進(jìn)行比對(duì)）
Python2 dexseq_prepare_annotation.py 基因組索引文件gtf? 扁平化文件.gtf
量化? python2 dexseq_count.py –p yes–s no –r name扁平化的文件.gtf 按照名字排序后的sam文件 輸出文件.txt
（7）組裝：
???? Cufflinks -P 8 –O保存文件夾/?比對(duì)后的文件bam/
De novo:
???????? Trinity.pl –seqType fq –JM 10G –left 1.fq–-right 2.fq –CPU 4
（8）差異表達(dá)分析：R語(yǔ)言作圖 ? ??
參照：http://www.360doc.com/content/18/0309/18/33459258_735717104.shtml ? ? ? ?http://www.360doc.com/content/17/0911/22/46164085_686360709.shtml ??http://www.biotrainee.com/thread-1084-1-1.html
setwd("c:/Users/du/Desktop/R/trails/")
##建立數(shù)據(jù)框? 轉(zhuǎn)錄組數(shù)據(jù)比對(duì)量化結(jié)果
condition1<-read.table("59counts.txt",sep = "\t",col.names=c("gene_id","condition1"))
condition2<-read.table("61counts.txt",sep = "\t",col.names = c("gene_id","condition2"))
rep1<-read.table("60counts.txt",sep = "\t",col.names = c("gene_id","x1"))
rep2<-read.table("62counts.txt",sep = "\t",col.names = c("gene_id","x2"))
raw_count<-merge(merge(condition1,condition2,by="gene_id"),merge(rep1,rep2,by="gene_id"))
Now if you have some more new interesting idea you will jump it > for example you can use join 1.txt 2.txt >1_2.txt ?Then you will get a mature txt and can be used in frame.OK let try it latter.
raw_count_filter<-raw_count[-1:-5,]
ENSEMBL <- gsub('\\.\\d*', '', raw_count_filter$gene_id)
row.names(raw_count_filter) <- ENSEMBL
raw_count_filter <- raw_count_filter[ ,-1]
##設(shè)計(jì)矩陣
group<-as.data.frame(c(rep("condition",2),rep("x",2)))
rownames(group)<-names(raw_count_filter)
names(group)<-"m"
#標(biāo)準(zhǔn)化矩陣
dds<-DESeqDataSetFromMatrix(countData = raw_count_filter,colData = group,design = ~m)
dds<-DESeq(dds)
summary(results(dds))
res<-results(dds,contrast<-c("m","condition","x"))
ressult<-res[order(res$padj),]
write.table(ressult,file="gene1_RNAseq2.xls",sep = "\t",quote=F,row.names = T,col.names = NA)
###提取差異分析結(jié)果
x<-read.table("gene1_RNAseq2.xls",sep = "\t",header = T,row.names = 1)
x
attach(x)
x1<-subset(x,padj<0.05 &(log2FoldChange>=1 | log2FoldChange <=-1))
write.table(x1,"differ_out_gene.xls",sep = "\t",quote=F,row.names = T,col.names = NA)
##GO KEGG GSEA analysis? ***（clusterProfiles package 可以下載現(xiàn)成的包，再導(dǎo)入熔任。)? 還需要DOSE? DO.db包
##小鼠的如下：? 注釋數(shù)據(jù)在R包中有
library(org.Mm.eg.db)
keytypes(org.Mm.eg.db)
library(clusterProfiler)
library(DOSE)
#GO analysis
#gene_id 轉(zhuǎn)換
gene<-row.names(x1)
transsid<-select(org.Mm.eg.db,keys=gene,columns=c("GENENAME","SYMBOL","ENTREZID"),keytype="ENSEMBL")
#看單個(gè)基因的話
single_Ap1m2<-select(org.Mm.eg.db,keys = "Ap1m2",columns = c("ENSEMBL","ENZYME","GENENAME", "SYMBOL","PROSITE","PMID","UNIGENE","UNIPROT","GO"),keytype = "SYMBOL")
##**正式GO分析
ego<-enrichGO(gene = row.names(x1),OrgDb = org.Mm.eg.db,keyType ="ENSEMBL",ont="MF" )
###ego1<-enrichGO(gene = row.names(x1),OrgDb = org.Mm.eg.db,keyType ="ENSEMBL",ont="ALL", qvalueCutoff = 0.01)
summary(ego)
head(summary(ego))
? #gorich 后做幾個(gè)圖：
##網(wǎng)絡(luò)圖
emapplot(ego)
#goplot(ego)
goplot(ego)
#Bar plot
barplot(ego, showCategory=20)
# 氣泡圖
dotplot(ego,font.size=5,showCategory=30)
#Gene-Concept Network
cnetplot(ego)
#UpSet Plot
upsetplot(ego)
#Heatmap-like functional classification
heatplot(ego)
###GSEA分析
#獲取按照l(shuí)og2FC大小來(lái)排序的基因列表
genelist <- x1$log2FoldChange
names(genelist) <- rownames(x1)
genelist <- sort(genelist, decreasing = TRUE)
# GSEA分析（具體參數(shù)參考：https://mp.weixin.qq.com/s/p-n5jq5Rx2TqDBStS2nzoQ）
gsemf<-gseGO(genelist,OrgDb = org.Mm.eg.db,keyType = "ENSEMBL",ont = "MF")
head(gsemf)
gseaplot(gsemf,geneSetID = "GO:0000977")
###KEGG（pathway）分析
# 轉(zhuǎn)換ID適合KEGG
m=bitr(rownames(x1),fromType = 'ENSEMBL',toType = 'ENTREZID', OrgDb = 'org.Mm.eg.db')
kegg <- m[,2]
# KEGG分析褒链，在KEGG官網(wǎng)中，物種都有對(duì)應(yīng)的縮寫(xiě)疑苔，小鼠mmu甫匹，其他的縮寫(xiě)看官網(wǎng)：http://www.genome.jp/kegg/catalog/org_list.html
kk <- enrichKEGG(kegg, keyType = 'kegg',organism = 'mmu', pvalueCutoff = 0.05, pAdjustMethod = 'BH', qvalueCutoff = 0.1)
head(summary(kk))
# 氣泡圖
dotplot(kk, font.size=5)
# 將GO/KEGG結(jié)果轉(zhuǎn)換成CSV格式輸出
write.table(as.data.frame(kk),"KEGG-enrich.xls",row.names =F)
write.csv(as.data.frame(ego),"GO-enrich.xls",row.names =F)
NOTICE：最近在學(xué)習(xí)和總結(jié)以及做一些自己的課題，因此會(huì)時(shí)常更新一部分惦费，敬請(qǐng)?jiān)彛?/p>
我看了一下：HTSeq軟件與cufflinks差異區(qū)別在哪里
如下：HTSeq產(chǎn)生的是reads匹配到基因外顯子上的序列兵迅，得到的是數(shù)據(jù)庫(kù)中基因名字（需要在R中轉(zhuǎn)換為俗稱(chēng)可認(rèn)識(shí)的名字，之后獲得差異表達(dá)基因）圖一薪贫；cufflinks：如果直接生成的是有FPKM值跟HTSeq類(lèi)似吧圖二恍箭；而cufflinks先組裝再合并再獲得差異基因，之后提取差異基因結(jié)果就是免去了R語(yǔ)言中很多的操作圖三與圖四瞧省；如下：
后續(xù)理解上：
事實(shí)上：可以用hisat2代替tophat2扯夭，畢竟作者也是這么建議的。
步驟分為如下過(guò)程臀突，為了好理解：
1. 數(shù)據(jù)下載：geo數(shù)據(jù)庫(kù)的測(cè)序數(shù)據(jù)勉抓；UCSC網(wǎng)站基因組數(shù)據(jù)（chromFa.tar.gz）；gencode網(wǎng)站基因組注釋文件（gtf）候学；hisat2網(wǎng)站的index文件藕筋；RSeQc軟件網(wǎng)站作覆蓋度的文件（.bed）
2.數(shù)據(jù)比對(duì)：hisat2 ?只是用了index文件； samtools的格式轉(zhuǎn)換--排序等
3.比對(duì)結(jié)果質(zhì)檢：RSeQc的質(zhì)檢梳码；
4.read計(jì)數(shù)隐圾，read歸類(lèi)于哪個(gè)基因（HTSeq）、轉(zhuǎn)錄（cufflinks）掰茶、外顯子區(qū)域（DEXSeq）暇藏；
5. 差異表達(dá)分析：歸類(lèi)后要進(jìn)行比較了；DESeq2包
6. 富集分析：因?yàn)槟骋粋€(gè)基因不能僅憑借表達(dá)量多少就判斷多少了濒蒋，如果是低表達(dá)的突然高表達(dá)一點(diǎn)呢盐碱，高表達(dá)但是比正常狀態(tài)卻少了呢把兔，這時(shí)候就需要看看富集到一起的結(jié)果是咋樣的了。Y叔的clusterfiler包瓮顽。GO KEGG
7. 其他分析：聚類(lèi)分析圖县好；主成分分析；（R語(yǔ)言中高級(jí)技能包括四類(lèi)：廣義線性模型暖混、聚類(lèi)分析缕贡、時(shí)間序列、主成分分析）拣播；別的嘛晾咪，需要什么就怎么操作吧。
