Bulk RNA-Seq轉錄組學習

與之對應的是single cell RNA-Seq，后面也會有類似文章。

參考：https://github.com/xuzhougeng/Learn-Bioinformatics/

作業：RNA-seq基礎入門傳送門

資料：RNA-seq Data Analysis-A Practical Approach(2015)

Bioinformatic Data Skill

biostar handbook

A survey of best practices for RNA-seq data analysis

Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis

Detecting differential usage of exons from RNA-seq data

轉錄組入門（1）: 工作環境准備

miniconda2

cd src
wget https://repo.continuum.io/miniconda/Miniconda2-latest-Linux-x86_64.sh
bash Miniconda2-latest-Linux-x86_64.sh

conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/msys2/
conda config --add channels bioconda
conda config --set show_channel_urls yes

在里面安裝各種工具

conda create -n biostar sra-tools fastqc hisat2 samtools htseq

 

轉錄組入門(2)：讀文章拿到測序數據

AKAP95 regulates splicing through scaffolding RNAs and RNA processing factors. Nat Commun 2016 Nov 8;7:13347. PMID: 27824034

GSE81916，https://www.ncbi.nlm.nih.gov/geo/

for i in `seq 48 62`;
do
    prefetch SRR35899${i}
done

用到RNA-seq的部分：

AKAP95 mainly promotes inclusion of exons globally.

To assess the global impact of AKAP95 on AS, we depleted 耗盡 AKAP95 in human 293 cells and mouse ES cells, and performed RNA-seq followed by DEXseq analysis29 to identify the differential exon usage in the cellular mRNAs.

（DEXSeq包是用來分析RNA-seq實驗數據中exon usage的差異，這里exon usage的差異指的是由於實驗條件導致的相對不同的exon usage）

文章的整體邏輯：

 

轉錄組入門(3)：了解fastq測序數據

for id in `seq 56 62`
do
    fastq-dump --gzip --split-3 -O /mnt/f/Data/RNA-Seq -A SRR35899${id}
done

FastQC質量報告

fastqc seqfile1 seqfile2 .. seqfileN
常用參數：
-o： 輸出路徑
--extract: 輸出文件是否需要自動解壓 默認是--noextract
-t: 線程， 和電腦配置有關，每個線程需要250MB的內存
-c: 測序中可能會有污染， 比如說混入其他物種
-a: 接頭
-q: 安靜模式

zcat SRR3589956_1.fastq.gz | fastqc -t 4 stdin
fastqc SRR3589956_1.fastq.gz

conda install multiqc
multiqc --help

轉錄組入門（4）：了解參考基因組及基因注釋

cd /mnt/f/Data
mkdir reference && cd reference
mkdir -p genome/hg19 && cd genome/hg19
nohup wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz &
tar -zvf chromFa.tar.gz
cat *.fa > hg19.fa
rm chr*

nohup wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_26/GRCh37_mapping/gencode.v26lift37.annotation.gtf.gz &
nohuop wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_26/GRCh37_mapping/gencode.v26lift37.annotation.gff3.gz &

 

轉錄組入門(5)： 序列比對

RNA-Seq數據分析分為很多種，比如說找差異表達基因或尋找新的可變剪切。如果找差異表達基因單純只需要確定不同的read計數就行的話，我們可以用bowtie, bwa這類比對工具，或者是salmon這類align-free工具，並且后者的速度更快。

但是如果你需要找到新的isoform，或者RNA的可變剪切，看看外顯子使用差異的話，你就需要TopHat, HISAT2或者是STAR這類工具用於找到剪切位點。因為RNA-Seq不同於DNA-Seq，DNA在轉錄成mRNA的時候會把內含子部分去掉。所以mRNA反轉的cDNA如果比對不到參考序列，會被分開，重新比對一次，判斷中間是否有內含子。

基本上就是HISAT2和STAR選一個就行。

cd referece && mkdir index && cd index
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz
tar -zxvf hg19.tar.gz

# 其實hisat2-buld在運行的時候也會自己尋找exons和splice_sites，但是先做的目的是為了提高運行效率
extract_exons.py gencode.v26lift37.annotation.sorted.gtf > hg19.exons.gtf &
extract_splice_sites.py gencode.v26lift37.annotation.gtf > hg19.splice_sites.gtf &
# 建立index， 必須選項是基因組所在文件路徑和輸出的前綴
hisat2-build --ss hg19.splice_sites.gtf --exon hg19.exons.gtf genome/hg19/hg19.fa hg19

mkdir -p RNA-Seq/aligned
for i in `seq 57 62`
do
    hisat2 -t -x reference/index/hg19/genome -1 RNA-Seq/SRR35899${i}_1.fastq.gz -2 SRR35899${i}_2.fastq.gz -S RNA-Seq/aligned/SRR35899${i}.sam &
done

# 運行實例

hisat2-2.0.4/hisat2 --no-discordant -t -x Database/hg19/GenomeHisat2Index/chrALL -U split_read.1.fq 2> Map2GenomeStat.xls | samtools view -bS - -o hisat2.bam

mkdir -p RNA-Seq/aligned
for i in `seq 56 58`
do
    hisat2 -t -x reference/index/hg19/genome -1 RNA-Seq/SRR35899${i}_1.fastq.gz -2 SRR35899${i}_2.fastq.gz -S RNA-Seq/SRR35899${i}.sam &
done

for i in `seq 56 58`
do
    samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam
    samtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bam
    samtools index SRR35899${i}_sorted.bam
done

比對質控(QC)

Picard https://broadinstitute.github.io/picard/
RSeQC http://rseqc.sourceforge.net/
Qualimap http://qualimap.bioinfo.cipf.es/

# Python2.7環境下
pip install RSeQC

# 先對bam文件進行統計分析， 從結果上看是符合70~90的比對率要求。

bam_stat.py -i SRR3589956_sorted.bam

# 基因組覆蓋率的QC需要提供bed文件，可以直接RSeQC的網站下載，或者可以用gtf轉換
read_distribution.py -i RNA-Seq/aligned/SRR3589956_sorted.bam -r reference/hg19_RefSeq.bed

 

轉錄組入門(6)： reads計數

如果你只需要知道已知基因的表達情況，那么可以選擇alignment-free工具（例如salmon, sailfish)，如果你需要找到noval isoforms，那么就需要alignment-based工具（如HISAT2, STAR）。

定量分為三個水平

基因水平(gene-level)
轉錄本水平(transcript-level)
外顯子使用水平(exon-usage-level)。

計數分為三個水平： gene-level, transcript-level, exon-usage-level

標准化方法： FPKM RPKM TMM TPM

# 安裝
conda install htseq
# 使用
# htseq-count [options] <alignment_file> <gtf_file>
htseq-count -r pos -f bam RNA-Seq/aligned/SRR3589957_sorted.bam reference/gencode.v26lift37.annotation.sorted.gtf > SRR3589957.count

mkdir -p RNA-Seq/matrix/
for i in `seq 56 58`
do
    htseq-count -s no -r pos -f bam RNA-Seq/aligned/SRR35899${i}_sorted.bam reference/gencode.v26lift37.annotation.sorted.gtf > RNA-Seq/matrix/SRR35899${i}.count 2> RNA-Seq/matrix/SRR35899${i}.log
done

paste *.txt | awk '{printf $1 "\t";for(i=2;i<=NF;i+=2) printf $i"\t";printf $i}'

 

轉錄組入門(7)： 差異表達分析

 

 

轉錄組入門(8)： 富集分析