Bulk RNA-Seq转录组学习
与之对应的是single cell RNA-Seq，后面也会有类似文章。
参考：https://github.com/xuzhougeng/Learn-Bioinformatics/
资料：RNA-seq Data Analysis-A Practical Approach(2015)
Bioinformatic Data Skill
A survey of best practices for RNA-seq data 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)： 差异表达分析
