RNA-seq 转录组数据分析处理（上）¶
1.背景知识补充¶
1.1 什么是RNA-seq？¶
RNA-seq即RNA-sequencing，又称转录组测序技术.它是一种通过观察基因表达来分析整个基因组的技术。主要测序对象包括信使RNA（mRNA）、微RNA（smallRNA）和非编码RNA（noncoding RNA）等或者其中一些，用高通量测序技术进行测序分析，反映出它们的表达水平.
1.2 高通量测序技术（high-throughput sequencing）¶
顾名思义，就是通量很高（对比sanger测序）的测序，又称“下一代”测序技术（Next-generation sequencing technology）（NGS），一次对几十万到几百万条DNA分子进行序列测定，同时高通量测序使得对一个物种的转录组和基因组进行细致全貌的分析成为可能，所以又被称为深度测序（deep sequencing).
1.3 RNA-seq的应用¶
RMA测序不但可以对基因表达的情况进行整理和归类，也可以对基因表达情况进行比较.在已有的研究中认为，一个物种个体的所有细胞里含有的DNA应该是一模一样的，但是每一种细胞里所表达的RNA都有细微的差异，这些细胞才会具有不同的功能。这也可以解释为什么肿瘤细胞会与机体细胞不同.通过RNA测序，可以了解许多基因的转录产物，进而可以更好的治疗疾病.
2.RNA-seq主要步骤¶
2.1 构建测序文库¶
分离RNA-->将RNA打断成小片段-->将小RNA片段反转录成DNA-->加接头-->PCR扩增（只有加上接头的片段才能被扩增）-->质量检查QC（看下文库的浓度和片段的长度）
注：我们建的文库用用的是cDNA而不是rna. 除了单链RNA不稳定外，还有一小部分原因是DNA的建库流程已经确定了，只要把RNA变成DNA后面的流程完全一样。
2.2 对测序文库进行测序¶
一块测序板上能包含多于400,000,000个片段，垂直于测序板排列。 测序仪有四种颜色的荧光探针A、T、C、G，与测序片段上碱基互补，结合上就闪一下自己带的荧光，比如A带红光，G带蓝光，C绿光，T橙光。许许多多的测序片段中同一排的碱基测完了，就把原来荧光的那个碱基冲掉了，再放下一个荧光碱基进来结合、放光，测序就是这样，结果就得到了raw data，就是fastq数据。
3.准备工作¶
学习illumina公司测序原理
测序的到的fastq文件
注释文件和基因组文件的准备
3.1 fastq文件的获取¶
步骤为：从数据库下载SRA数据-->SRA转化为fastq文件
数据库资源： NCBI数据库： https://www.ncbi.nlm.nih.gov/
打开GEO数据库，输入GSE号.Download family是经过GEO数据库整理和标准化的数据.Supplementary file是原始数据，提供SRA格式的文件.
打开SRA数据库，输入SRR号
3.2 简单介绍fastq文件¶
在illumina测序文件中，采用双端测序（paired-end），一个样本得到的是seq_1.fastq.gz和seq_2.fastq.gz,每个文件存放一段测序文件.
在illumina公司测得的序列文件经过处理以fastq文件协议储存为.*fastq格式文件.在fastq文件中每4行储存一个read.
第一行：第1行主要储存序列测序时的坐标等信息.
第二行：第2行就是测序得到的序列信息，一般用ATCGN来表示，其中N表示荧光信号干扰无法判断到底是哪个碱基.
第三行：第3行以“+”开始，可以储存一些附加信息，一般是空的.
第四行：第4行储存的是每个碱基的质量得分，与第2行的碱基序列是一一对应的，其中的每一个符号对应的ASCII值成为phred值，可以简单理解为对应位置碱基的质量值，越大说明测序的质量越好。用ASCII码表示碱基质量是为了减少文件空间占据和防止移码导致的数据损失.
如果要仔细了解fastq文件，请打开链接https://zhuanlan.zhihu.com/p/20714540
3.2 注释文件（gtf/gff）和基因组文件(fasta)的准备¶
注释文件下载：
wet ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M25/gencode.vM25.chr_patch_hapl_scaff.annotation.gtf.gz
基因组文件下载：参考基因组下载：服务器上已经下载好 /public/home/zhengxq/genome，无需下载.hg38相当于hg19的一个重大的进步
4.RNA-seq处理步骤¶
1.SRA数据下载：prefetch
2.SRA转为fastq格式: fastq-dump
3.质控: fastqc,multiqc
4.去杂躁音（去除低质量的adapter）: trim_galore
5.基因组比对:hisat2
6.sam文件转bam文件:samtools
7.对bam文件按按read name排序：samtools
8.reads计数,合并矩阵并进行注释：htseq-count
4.1 SRA数据下载¶
我们对4个来自小鼠的测序数据进行分析：
prefetch SRR3589959
prefetch SRR3589960
prefetch SRR3589961
prefetch SRR3589962
注：下载完成后，文件存储在：/public/home/zhengxq/ncbi/public/sra，可复制到自己的文件夹下.
4.2 SRA转为fastq格式: fastq¶
for ((i=59;i<=62;i++));do fastq-dump --gzip --split-3 -A SRR31621$i.sra -O .;done
4.3.检查下载的数据是否完整：md5值¶
md5sum *gz >md5.txt
cat md5.txt
md5sum -c md5.txt
4.4 质控: multiqc¶
fastqc -o . *.fastq.gz
multiqc .
5.基因组比对:hisat2¶¶
建立hisat2索引
nohup hisat2-build -p 4 mm10.fa genome &
命令完成后，会在mm10文件夹下生成如下文件：
开始比对：用hisat2，得到SAM文件¶
我的fastq文件在/public/home/zhengxq/zhuyj/RNAseq2/rawdata
我的reference在/public/home/zhengxq/zhuyj/RNAseq2/rawdata/mm10
我的index在/public/home/zhengxq/zhuyj/RNAseq2/rawdata/mm10
比对后得到的bam文件会存放在/public/home/zhengxq/zhuyj/RNAseq2/rawdata/aligned
source ~/miniconda3/bin/activate
for ((i=59;i<=62;i++));
do hisat2 -t -x /public/home/zhengxq/zhuyj/RNAseq2/rawdata/mm10/genome -1 /public/home/zhengxq/zhuyj/RNAseq2/rawdata/SRR35899${i}.sra_1.fastq.gz -2 /public/home/zhengxq/zhuyj/RNAseq2/rawdata/SRR35899${i}.sra_2.fastq.gz -S SRR35899${i}.sam ;
done
6. SAM文件转换为BAM文件¶
cd /public/home/zhengxq/zhuyj/RNAseq2/rawdata/aligned
for ((i=59;i<=62;i++));do samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam;done
7. bam文件按reads name排序¶
for ((i=59;i<=62;i++));
do samtools sort -n SRR35899${i}.bam SRR35899${i}_nsorted.bam;
done
8.reads计数，合并矩阵并进行注释¶
注：先解压注释文件，注释文件为gencode.vM25.chr_patch_hapl_scaff.annotation.gtf.gz
for ((i=59;i<=62;i++));
do htseq-count -r name -f bam /public/home/zhengxq/zhuyj/RNAseq2/rawdata/aligned/SRR35899${i}_nsorted.bam /public/home/zhengxq/zhuyj/RNAseq2/rawdata/mm10/gencode.vM25.chr_patch_hapl_scaff.annotation.gtf > /public/home/zhengxq/zhuyj/RNAseq2/rawdata/count/SRR35899${i}.count;
done
wc -l *.count
看下每个文件的格式，查看前4行,第一列ensembl_gene_id,第二列read_count计数¶
head -n 4 *.count
后4行¶
tail -n 4 *.count
8.1 合并表达矩阵并进行基因名注释（Rstudio中进行）¶
Tumor:SRR316214,SRR316215
Adjacent Normal Liver:SRR316212,SRR316213
options(stringsAsFactors = FALSE)
control1<-read.table("SRR3589959.count",sep = "\t",col.names = c("gene_id","control1"))
head(control1)
control2<-read.table("SRR3589961.count",sep = "\t",col.names = c("gene_id","control2"))
treat1<-read.table("SRR3589960.count",sep = "\t",col.names = c("gene_id","treat1"))
treat2<-read.table("SRR3589962.count",sep = "\t",col.names = c("gene_id","treat2"))
数据整合merge¶
gencode的注释文件中的gene_id（如ENSMUSG00000105298.13_3）在EBI是不能搜索到的，所以用gsub功能只保留ENSMUSG00000105298这部分
处理之前先看一下,也就是最后5行是我们不需要的，可以删除¶
tail(control1)
进行整合¶
raw_count <- merge(merge(control1, control2, by="gene_id"), merge(treat1, treat2, by="gene_id"))
head(raw_count)
tail(raw_count)
删除前5行¶
raw_count_filt <- raw_count[-1:-5,]
head(raw_count_filt)
因为我们无法在EBI数据库上直接搜索找到ENSMUSG00000024045.5这样的基因，只能是ENSMUSG00000024045的整数，没有小数点，所以需要进一步替换为整数的形式。
# 第一步将匹配到的.以及后面的数字连续匹配并替换为空，并赋值给ENSEMBL
ENSEMBL <- gsub("\\.\\d*", "", raw_count_filt$gene_id)
# 将ENSEMBL重新添加到raw_count_filt1矩阵
row.names(raw_count_filt) <- ENSEMBL
head(raw_count_filt)
保存count矩阵，方便后续进行下游分析¶
write.csv(raw_count_filt, file='raw_count_filt.csv')