|
本文是由比利時(shí)列日大學(xué)Marc HANIKEN課程整理,。陸陸續(xù)續(xù)交付NE大學(xué)??1個(gè)月完成，根據(jù)需要做的內(nèi)容分為四個(gè)部分,。 1 目標(biāo)RNA-Seq的目標(biāo)是說(shuō)明如何處理和分析RNA-Seq的數(shù)據(jù)以識(shí)別差異基因（DGE）,。 2 數(shù)據(jù)介紹擬植物南芥的基因型（wt 組織模式和多種體型）在（c）和處理（t）條件下處于下,。本樣品獨(dú)立株植物實(shí)驗(yàn)） 3 3 次 NextSeq 0 儀器以 4 次重復(fù)使用Illumina,。集群基因組和芯片組在整個(gè)組件中運(yùn)行，使用中和兩個(gè) bp 的 5 個(gè)端快速介紹,。 ,。 3 分析數(shù)據(jù)3.1 查看數(shù)據(jù)head <your_sample>.fastq  圖像.png
將所有的樣本名稱攔截一個(gè)文件 for f in *.fastq; do echo `basename $f .fastq`; done > samples.ids
3.2 RNA-Seq 數(shù)據(jù)分析中讀取映射的一般考慮
第一部分：3.3 讀取映射到參考基因組3.3.1 工具介紹1. 頂帽軟件我們將使用流行的帽子，這是將 RNA-Seq 外接閱讀與基因組外顯子以識(shí)別子程序-顯式程序的短連接的,。更多 TopHat 如何找到連接點(diǎn)的原理： TopHat 可以通過(guò)注釋的情況下將 RNA-Seq 讀取到?jīng)]有參考基因。這個(gè)映射信息,，TopHat建立一個(gè)可能的剪接連接的數(shù)據(jù)庫(kù),，然后將讀取映射到這些連接以確認(rèn)它們。 這一段讀到這個(gè)標(biāo)題的機(jī)子可能有1個(gè)00個(gè)遺漏或短問(wèn)題的外顯,，但會(huì)在最初的象征中將比所有的內(nèi)容都被更多地漏掉,。獨(dú)立映射這些。 TopHat 兩個(gè)約定生成可能的剪接點(diǎn)數(shù)據(jù)庫(kù),。這種情況,，“GTAG”、“GC-AG”和“AT-AG”和“AT-AG”和“AT-AG” AC“通常在其中含有不同品種的標(biāo)題尋找,。第二個(gè)來(lái)源是“封面年齡的島嶼”的開(kāi)始,，是最終中部的中部地區(qū)尋找到的。將這些內(nèi)含子連接起來(lái)的方法,。我們只建議第二個(gè)選項(xiàng)（--coverage-search）用于將短讀?。?lt;45bp）和用戶讀?。?lt;=1000萬(wàn)）。后一個(gè)選項(xiàng)對(duì)“GT-AG”內(nèi)含子之間的比,。 Tophat可以使用FASTA,FASTQ(推薦)格式的讀取,。 想要使用這個(gè)軟件，首先需要使用一下命令：  圖像.png
Bowtie2用于熱門組上的閱讀,。 蝶領(lǐng)結(jié)擅長(zhǎng)使用一種超高配的技術(shù),，用于與組合工具和排列組合。 Bowtie 2 保持珠寶形狀使用組合（基于Browtie 2 對(duì)BWT 進(jìn)行）,，通常其占用或占用的內(nèi)存大小,。 Bowtie 2 的結(jié)構(gòu)需要占用多少個(gè)內(nèi)存。雙端模式,。同時(shí)可以使用多個(gè)處理器來(lái)更高的關(guān)注度,。 Bowtie 2 以SAM 格式輸出的其他方式，以SAM格式輸出的其他方式,，使用授權(quán)文件和大量使用同樣的工具（SAMtools,、GATK 的許可互操作）。Bowtie 2GPLv3 在和下分發(fā),，Mac OS X Linux BSD 和它在Windows 下的運(yùn)行,。 Bowtie 2和Bowtie 2和Bowtie BS （這里也叫“集成1 ” sowtie 2和Bowtie BS）通常是比較多種其他學(xué)組的，包括變異,、RNA-seq,、Ch IPeq。工具中,，這里有其中一些,。 要與 Tophat 的連接點(diǎn)，您首先需要為 RNA-Seq 中的生物體安裝蝴蝶結(jié)指數(shù),。使用 bowtie2-build 很容易自己制造一個(gè),。  圖像.png
Bowtie2 從 bowtie 索引中提取信息，允許確定它是什么索引以及使用什么序列來(lái)制造它,。 2. GFF/GTF 格式文件通過(guò)基因特征（例如外含子/內(nèi)含子描述格式組的基因組）提供的基因組注釋文件,，可以幫助通過(guò)頂帽在基因組上進(jìn)行讀取映射。 注釋文件以 GFF/GTF 提供,。 Tophat 使用的基因組注釋文件就是 GFF/GTF  圖像.png
GTF(general transfer format)是GFF第二個(gè)版本， 3 htseq-count軟件給定一個(gè)具有組合范圍的基因的文件,，htseq-count 會(huì)計(jì)算出有多少讀取的特征映射到某個(gè)特征列表,。 - 在每個(gè)情況下，特征通常是每個(gè)基因被結(jié)合的,，其中所有外顯子的地方也可以顯示子的一個(gè)特征,，例如,，為了檢查。對(duì)于比較 ChIP-Seq,，特征可能是列表中的結(jié)合區(qū)域,。 htseq-count 腳本允許在不同模式之間進(jìn)行選擇。 hts-count 的位置重疊模式的工作原理如下：定義一個(gè)集合 S(i) 的位置為我重疊的特征的集合,。然后，考慮集合 S,，它是（我遍歷或讀取對(duì)中的所有位置）
3.3.2 下載擬南芥參考網(wǎng)址：https://www./（需要注冊(cè)） 3.3.3 給基因參考建索引使用bowtie2-build,。為擬南芥編制索引,，花費(fèi)2分鐘 bowtie2-build Arabidopsis.fasta At_ref
檢查指數(shù)，幾秒鐘 3.3.4 讀取映射內(nèi)容為存在以逗號(hào)隔FA打開(kāi)的FASTQ或STA格式文件 使用tophat完成一般使用命令：  圖像.png
更多的選擇閱讀文檔 其中： --num-threads 4 ##可以多線程 請(qǐng)注意,，所提供的 GTF/GFF 文件的第一個(gè)索引（指示特征所在的染色體相列或重疊群的列）的必須與 TopH 的 Bowtie-中的參考值,。您可以使用序列匹配檢查進(jìn)行 開(kāi)始操作軟參考鏈接組基因的FASTA： ln -s Arabidopsis.fasta At_ref.fa
創(chuàng)建簡(jiǎn)單的索引,。立即創(chuàng)建，方便使用所有樣本,，簡(jiǎn)單組圖 5 分鐘 會(huì)在transcriptome_data/下產(chǎn)生10個(gè)文件 映射閱讀,，先創(chuàng)建一個(gè)模板 tophat -o output_[% basename %] --read-mismatches 2 --min-intron-length 40 --max-intron-length 2000 --num-threads 2 --report-secondary-alignments --no-novel-juncs --transcriptome-index=transcriptome_data/At_ref At_ref [% basename %].fastq
樣品創(chuàng)建一個(gè)灰 提交任務(wù)： for f in `cat samples.ids`
do qsub -pe snode 2 tophat_$f.sh
done
此步驟費(fèi)用大約 1 小時(shí) 對(duì)所有的樣本進(jìn)行總結(jié)查看 for f in `cat samples.ids`
do head output_$f/align_summary.txt
done
3.3.5 讀計(jì)數(shù)使用htseq-count 圖像.png
指定輸出任務(wù)的一個(gè)表，包含功能（這里是由于計(jì)算）的計(jì)數(shù),，然后是特定測(cè)點(diǎn)的特殊點(diǎn),，用于未針對(duì)特定原因進(jìn)行的讀取。于過(guò)濾,。情況是：  圖像.png
提示：如果你有特定于鏈的特定數(shù)據(jù),，否則請(qǐng)確保你設(shè)置的 RNA-Seq 數(shù)據(jù)不是特定鏈的協(xié)議。-strand=no,！ htseq-count 有很多選項(xiàng),，請(qǐng)查看鏈接文檔 的 一些選項(xiàng)： -f < sam or bam># 輸入文件，sam 或 bam 格式 -s <yes/no/reverse> 讀計(jì)數(shù)模板運(yùn)行花費(fèi)半個(gè)小時(shí),。 搜索征集統(tǒng)計(jì)信息貝殼命令 for f in <your_name>_htseqcount_*.o*; do tail -n 5 $f; done
.組件計(jì)算矩陣基因的名字 識(shí)數(shù) for f in `cat samples.ids`
do cut -f2 <your_name>_htseqcount_$f.o* > $f.count
done
組件列表和計(jì)數(shù)得到這個(gè)結(jié)果文件,，將用于 GE 的統(tǒng)計(jì)分析,，第二部分： 4閱讀到參考組,。3.4.1 工具介紹
2組組您需要分析
分析使用每一個(gè)腳本：使用對(duì)齊的工具進(jìn)行統(tǒng)計(jì)分析,。因此，我們將使用對(duì)齊的工具來(lái)展示
3 3.4 擬擬南芥參考組2,。來(lái)自Araport，需要登錄進(jìn)行免費(fèi)注冊(cè),。再使用以下代碼獲取,。 curl -sO -H 'Authorization: Bearer <your_id_key>' https://api./files/v2/media/system/araport-public-files// \
Araport11_Release_201606/annotation/Araport11_genes.201606.cds.fasta.gz
3.4.索引擬南芥參考組3使用ltrinity的perl命令：align_and_estimate_abundance.pl，可以對(duì)所有樣本一次完成,。  圖像.png
索引的操作命令這個(gè)過(guò)程花費(fèi)大約5分鐘,，會(huì)生成14個(gè)文件，包含.bowtie2 .和.RSEM 3.4.4 對(duì)排列和計(jì)數(shù)使用 ltrinity 的 perl 命令：align_and_estimate_abundance.pl,，并使用 RSEM 估計(jì)方法  圖像.png
 圖像.png
2建立gene_trans_地圖 grep \> Arabidopsis_transcripts.fasta | cut -f2 -d '>' | cut -f1 -d '|' > transcripts.ids
# Let's paste twice this list in the same file
$ paste transcripts.ids transcripts.ids > double_transcripts.ids
$ head double_transcripts.ids
# And apply the following perl one liner to remove the transcript number
# from 1st column
$ perl -nle 's/^(AT\w+)\.\d+/$1/g; print' double_transcripts.ids > gene_trans_map.txt
3,、進(jìn)行地圖和計(jì)數(shù) align_and_estimate_abundance.pl 命令使用模板： 創(chuàng)建多個(gè)樣本的sh文件： for f in `cat samples.ids`
do tpage --define queue=smallnodes --define basename=$f --define thread=2 trinity_align_estimate.tt > align_estimate_$f.sh
done
提交任務(wù)： 這大概要花90分鐘 圖像.png
圖像.png
3.4.5 生成表達(dá)矩陣使用：trinity下的abundance_estimates_to_matrix.pl命令
將 perl /media/vol1/apps/trinityrnaseq-2.2.0/util/abundance_estimates_to_matrix.pl --est_method RSEM trinity_*/*.genes.results --out_prefix <your_name>
大概需要2分鐘 該腳本輸出多個(gè)文件 第三部分： 3.5 差異表達(dá)的基因使用R包DESeq2,。3.5.1 包介紹詳細(xì)文檔介紹：https:///packages/release/bioc/html/DESeq2.html,。 圖像.png
DESeq2將首先對(duì)數(shù)據(jù)進(jìn)行建模的例子,。 有關(guān)每個(gè)步驟的詳細(xì)信息，請(qǐng)參閱相應(yīng)手冊(cè)頁(yè)。調(diào)整值的信息,，請(qǐng)參見(jiàn)結(jié)果手冊(cè)頁(yè),。 使用DESeq(object)，是一個(gè)DESeqDataSet的對(duì)象,。如：DESeqDataSetFromMatrix,。 3 結(jié)果名稱 返回模型的估計(jì)模型（因子）的名稱
參數(shù)是DESeqDataSet已經(jīng)在其上調(diào)用中以下函數(shù)： DESeq 、bino值對(duì)比WaldTest或nbinomLRT之一,，對(duì)比值比較變化 alpha 優(yōu)化的顯著性結(jié)束值（默認(rèn)為 0.）。如果調(diào)整的 p 最終值 (FDR) 為 1,，則 alpha 應(yīng)設(shè)置為該值,。 3.5.2 下載DESeq2library(BiocManager)
BiocManager::install('openssl')
BiocManager::install('RCurl')
BiocManager::install(c('DESeq2','limma','gplots'), force = T)
3.5.3 特征基因表達(dá)差異（成對(duì)比較）我們將在下面發(fā)現(xiàn)的基因需要允許需要的 R 腳本,。您在里面按順序添加每個(gè)新步驟,。然后，根據(jù) DGE 的治療類型（Ctrl vs Treat）,，最后治療對(duì)各個(gè)種的類型,。基因中必須考慮到這一點(diǎn),。 Step 1. 加載數(shù)據(jù)并描述數(shù)據(jù)集 步驟 2. 建立基因型響應(yīng)分析模型 #Genotype effect
#####
#Load data using the DESeqDataSetFromMatrix command
genotDesign=DESeqDataSetFromMatrix(countData = countData,colData = colData,
design = ~ genot)
#Build model using the DESeq command
genot_DESeq <- DESeq(genotDesign)
#Observe parameters of the model
resultsNames(genot_DESeq)
步驟 3. 使用 PCA 對(duì)數(shù)據(jù)進(jìn)行匯總統(tǒng)計(jì) Step 4. 建立樣本距離的熱圖 #Build sample distance
sampleDist <- dist(t(assay(rld)))
#Build heatmap
sampleDistMatrix<-as.matrix(sampleDist)
rownames(sampleDistMatrix)<-paste(rld$g_t)
colnames(sampleDistMatrix)<-NULL
colours=colorRampPalette(rev(brewer.pal(9, 'Blues')))(300)
tiff(filename = 'heatmap_sampledist_Treat_root.tiff', width = 1500,
height = 1500, units = 'px', res = 150)
heatmap.2(sampleDistMatrix, dendrogram = 'both', trace = 'none', col = colours,
main = 'Treat Root Sample Distance', margin=c(6, 8))
dev.off()
步驟 5. 識(shí)別基因型主動(dòng)的 DGE 步驟 6,。 plotCounts(genot_DESeq, 'AT2G19110', intgroup = 'genot')
第四部分：3.6數(shù)據(jù)挖掘我們非常容易和我們一起使用 GE 數(shù)據(jù)集進(jìn)行的數(shù)據(jù)接口。 Thalemine 非常容易獲得相關(guān)數(shù)據(jù)集的功能,。 為了使用這個(gè),，我們首先需要從DESeq中 |
|
|
