轉錄組差異表達分析小實戰（一）
轉錄組差異表達分析小實戰（一）
讀文獻獲取資料
文獻名稱：AKAP95 regulates splicing through scaffolding
RNAs and RNA processing factors
-
查詢資料：Data availability
The RIP-seq an RNA-seq data have been deposited in the Gene
Expression Omnibus database, with accession code GSE81916. All other data is
available from the author upon reasonable request. -
獲得GSE號：GSE81916
下載測序資料
-
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE81916獲取資料資訊，並點選網址下方的ftp，下載測序資料
-
從https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA323422可知我們需要的mRNA測序編號為SRR3589956到SRR3589962
-
通過Apera下載SRR資料，這裡以SRR3589956為例：
ascp -T -i /home/anlan/.aspera/connect/etc/asperaweb_id_dsa.openssh[email protected]-private.ncbi.nlm.nih.gov:sra/sra-instant/reads/ByRun/sra/SRR/SRR358/SRR3589956/SRR3589956.sra ./
轉化fastq測序資料
-
通過sratoolkit工具將SRR檔案轉化為fastq格式的測序資料（寫了個shell迴圈）
for i in $(seq 56 62);do nohup fastq-dump --split-3 SRR35899${i} &;done -
通過fastqc對每個fastq檔案進行質檢，用multiqc檢視整體質檢報告（對當前目錄下的fastq測序結果進行質檢，生成每個fq檔案的質檢報告總multiqc整合後統計檢視）
fastqc *.fastq multiqc ./點選這個url可以檢視我這個multiqc報告：http://www.bioinfo-scrounger.com/data/multiqc_report.html
- 如果有接頭或者質量值不達標的需要進行過濾，這次的資料質量都不錯，因此直接進行比對即可
序列比對
-
安裝hisat2軟體，下載人類的hiast2索引檔案
-
hisat2下載並安裝：
ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat2-2.1.0-Linux_x86_64.zip unzip hisat2-2.1.0-Linux_x86_64.zip -
下載hisat2的human索引
ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz tar zxvf hg19.tar.gz
-
-
用hisat2進行比對，測序資料放在data目錄下，索引檔案放在reference/index/hisat2/hg19目錄下，SRR3589956-SRR3589958為人的測序資料
for i in $(seq 56 58);do hisat2 -p 4 \ -x ~/reference/index/hisat2/hg19/genome \ -1 ./data/SRR35899${i}_1.fastq -2 ./data/SRR35899${i}_2.fastq \ -S SRR35899$i.sam >SRR35899${i}.log;done -
用samtools將sam檔案轉化為bam檔案，並使用預設排序
for i in $(seq 56 58);do samtools sort [email protected] 5 -o SRR35899${i}.bam SRR35899${i}.sam;done
reads計數
-
用htseq對比對產生的bam進行count計數
-
htseq安裝，使用miniconda，省事！唯一的問題是htseq版本不是最新的，是0.7.2。想要最新版還是要正常安裝，可參考http://www.biotrainee.com/thread-1847-1-2.html
conda install -c bioconda htseq -
用htseq將對比後的結果進行計數
for i in $(seq 56 58);do htseq-count -f bam -r pos -s no \ SRR35899${i}.bam ~/reference/genome/hg19/gencode.v26lift37.annotation.gtf \ 1>SRR35899${i}.count 2>SRR35899${i}_htseq.log;done
-
-
將3個count檔案（SRR3589956.count，SRR3589957.count，SRR3589958.count）合併成一個count矩陣，這是就需要指令碼來解決這個問題，不然其他方法會稍微麻煩點
#!/usr/bin/perl -w use strict; my $path = shift @ARGV; opendir DIR, $path or die; my @dir = readdir DIR; my $header; my @sample; my %hash; foreach my $file (@dir) { if ($file =~ /^\w+.*\.count/) { push @sample, $file; $header .= "\t$file"; open my $fh, $file or die; while (<$fh>) { chomp; next if ($_ =~ /^\W+/); my @array = split /\t/, $_; $hash{$array[0]} -> {$file} = $array[1]; } close $fh; } } print "$header\n"; map{ my $gene = $_; print "$gene"; foreach my $file (@sample) { print "\t".$hash{$gene} -> {$file}; } print "\n"; }keys %hash; -
按照接下來的劇本，應該講count_matrix檔案匯入DESeq進行差異表達分析。但是從這篇文章的Bioinformatic analyses部分可以發現，作者的control組的2組資料是來自2個不同的批次（一個是SRR3589956，另外一個來源GSM1095127 in GSE44976），treat組倒是同一個批次（SRR3589957和SRR3589958）。但是對於Mouse cells來說，倒是滿足2個control和2個treat都正常來自同個批次，因此打算重新用SRR3589959-SRR3589962重新做個一個count_matrix進行後續差異分析