http://www.fungenomics.com/article/30 【專題】基因組學技術專題（二）—— 為什麼說FPKM/RPKM是錯的 下載資料 wget是linux下一個從網路上自動下載檔案的常用自由工具。它支援HTTP，HTTPS和FTP協議，可以使用HTTP代理。一般的使用方法是: wget + 空格 + 引數 + 要下載檔案的url路徑，例如： 1wget http://www.linuxsense.org/xxxx/xxx.tar.gz Wget常用引數 -b：後臺下載，Wget預設的是把檔案下載到當前目錄。 -O：將檔案下載到指定的目錄中。 -P：儲存檔案之前先建立指定名稱的目錄。 -t：嘗試連線次數，當Wget無法與伺服器建立連線時，嘗試連線多少次。 -c：斷點續傳，如果下載中斷，那麼連線恢復時會從上次斷點開始下載。 -r：使用遞迴下載 格式轉換 // use sra-tools to transform > fastq-dump *.sra 為了後面分析方便，把相應的序列檔名改成相應的組 mv SRR1871481.fastq WT_Rep1.fastq 質量控制 Pre-Alignment QC 使用fastqc 軟體來對原始序列fastq 檔案進行質量檢測 export PATH=/home/maque/FastQC/:$PATH fastqc *.fastq 這樣每個 fastq 檔案都能生成一個 html 報告檔案，很詳細 序列比對 使用 tophat 和 bowtie 進行比對 待新增 tophat2 -p 8 --bowtie1 -G genes.gtf -o WT_Rep1_output ../genome WT_Rep1.fastq 其他5個 fastq 檔案與上面一致 -p 8 使用8執行緒 --bowtie1 使用bowtie1 , 預設是bowtie2 -G 後面跟註釋檔案 gtf -o 後跟輸出資料夾 最後面跟 原始序列 fastq 檔案 經完成比對，生成了一個 accepted_hits.bam 檔案， 這個就是 samtools 生成的，後面主要也是利用這個檔案繼續分析。 建議提前利用 IGV 軟體檢視一下 該 bam 檔案，可以知道mapping 的情況。 檢視bam檔案 如果想檢視，需要先對 該bam檔案進行 index ,比如： samtools index WT_Rep1_output/accepted_hits.bam Use Cufflinks to generate expression estimates from the SAM/BAM files Use Cufflinks to generate expression estimates from the SAM/BAM files ref export PATH=/home/maque/cufflinks-2.2.1.Linux_x86_64/:$PATH cufflinks -p 8 -o WT_Rep1_cuffout WT_Rep1_output/accepted_hits.bam 其他5個與之類似 -p 8 使用8執行緒 -o WT_Rep1_cuffout 輸出目錄 最後面跟相應的 bam 檔案 該過程完成後，會生成相應的資料夾，裡面有相應的檔案，後面會使用 transcripts.gtf 檔案 Differential Expression ref ls -1 *cuffout/transcripts.gtf > assembly_GTF_list.txt cuffmerge -p 8 -o merged -g Arabidopsis_thaliana_Ensembl_TAIR10/Arabidopsis_thaliana/Ensembl/TAIR10/Annotation/Genes/genes.gtf -s Arabidopsis_thaliana_Ensembl_TAIR10/Arabidopsis_thaliana/Ensembl/TAIR10/Sequence/WholeGenomeFasta/genome.fa assembly_GTF_list.txt -p 8 使用8執行緒 -o merged 後跟目錄 -g 後跟參考基因的gtf檔案 -s 後跟基因組序列檔案 最後跟上一步新建的 assembly_GTF_list.txt 接下來使用 cuffdiff cuffdiff -o rna_de/diff_out -p 8 -L WT,athb merged/merged.gtf WT_Rep1_output/accepted_hits.bam,WT_Rep2_output/accepted_hits.bam,WT_Rep3_output/accepted_hits.bam athb_Rep1_output/accepted_hits.bam,athb_Rep2_output/accepted_hits.bam,athb_Rep3_output/accepted_hits.bam -o 後跟輸出檔案目錄 -p 8 使用8執行緒 -L WT,athb '-L' tells cuffdiff the labels to use for samples 後跟 上一步由 cuffmerge 生成的 merged.gtf 檔案 最後跟6個bam 檔案， 組內由逗號隔開，組間由空格隔開。 該過程會新建一個diff_out 資料夾，裡面包含了很多資訊 這些資訊可以使用 R 包 cummeRbund 很方便的進行分析 使用cummeRbund 文件 推薦流程 可以按照推薦流程文件中的步驟去分析 下面主要說一些注意點： 安裝 source("http://bioconductor.org/biocLite.R") biocLite("cummeRbund") 讀入資料 首先先 cd 到上一個步驟生成的 diff_out 資料夾 library(cummeRbund) cuff <- readCufflinks() 這樣即完成讀入資料。 各種分析圖 可以按照推薦流程中的去分析，也可以參考 Nature Protocol 尋找差異表達基因 大部分進行RNA-Seq 實驗的目的主要還是尋找實驗組與對照組之間的差異表達基因。 一種方法是: mySigGeneIds<-getSig(cuff,alpha=0.05,level='genes') mySigGeneIds length(mySigGeneIds) mySigGenes<-getGenes(cuff,mySigGeneIds) mySigGenes diffData(mySigGenes) featureNames(mySigGenes) 另一種方法是： gene_diff_data <- diffData(genes(cuff)) sig_gene_data <- subset(gene_diff_data, (significant == 'yes')) sig_gene_data 這些方法列出的 gene_id 是像 XLOC_000268 這樣的格式， 它所對應的通用的gene_id 比如AT1G06080 , 這種一一對應關係檔案可以在合併的 merged.gtf 檔案中尋找 而AT1G06080 這種gene_id 所對應的基因型別，基因名稱等資訊可以在 tair10 資料夾中的 gene.gtf 檔案中找到 AT1G06080 這種gene_id 所對應的基因名稱也可以在 同一資料夾中的 refFlat.txt 檔案中找到。 也可以先把上一步輸出的 gene_id 放到一個 list.txt 中，注意不要有空行，最好使用 vim , 然後利用 grep 即可： grep -f list.txt merged.gtf | less - S 以上就是rna-seq 資料分析的簡單過程，很多細節沒有提，而且還有很多其他步驟可以進行擴充套件，這些還需要再進一步深入理解。