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(看下文库的浓度和片段的长度)

3b7cce8d700191a56a5b01205ac5a99.png

注:我们建的文库用用的是cDNA而不是rna. 除了单链RNA不稳定外,还有一小部分原因是DNA的建库流程已经确定了,只要把RNA变成DNA后面的流程完全一样。

2.2 对测序文库进行测序

一块测序板上能包含多于400,000,000个片段,垂直于测序板排列。 测序仪有四种颜色的荧光探针A、T、C、G,与测序片段上碱基互补,结合上就闪一下自己带的荧光,比如A带红光,G带蓝光,C绿光,T橙光。许许多多的测序片段中同一排的碱基测完了,就把原来荧光的那个碱基冲掉了,再放下一个荧光碱基进来结合、放光,测序就是这样,结果就得到了raw data,就是fastq数据。

b5c4b149840abf385be129e91642bc9.png

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码表示碱基质量是为了减少文件空间占据和防止移码导致的数据损失.

228e20da8f41e91d30f2fa6f5b0950a.png

如果要仔细了解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值

In [ ]:
md5sum *gz >md5.txt
cat md5.txt

01a348af59dcbe252470631d1ba7001.png

In [ ]:
md5sum -c md5.txt

b6f405587673d1de8739b3714f24f6c.png

4.4 质控: multiqc

In [ ]:
fastqc -o .  *.fastq.gz
In [ ]:
multiqc .

5.基因组比对:hisat2¶

建立hisat2索引

nohup hisat2-build -p 4 mm10.fa genome &

命令完成后,会在mm10文件夹下生成如下文件:

7acbde92aca79f7f546aa3ee897391a.png

开始比对:用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

In [ ]:
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

81192f99cc6c8c24b53545a1f6b84b6.png

6. SAM文件转换为BAM文件

In [ ]:
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

cdb68ea8ffefad2958e94c2b1dba86d.png

7. bam文件按reads name排序

In [ ]:
for ((i=59;i<=62;i++));
do samtools sort -n SRR35899${i}.bam SRR35899${i}_nsorted.bam;
done

61f53b4ff3a7e26641e9a5667d334b1.png

1e343d44aa2957d5c573b389f8f392f.png

8.reads计数,合并矩阵并进行注释

注:先解压注释文件,注释文件为gencode.vM25.chr_patch_hapl_scaff.annotation.gtf.gz

In [ ]:
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

45dd99be4175e575d2b19cdd926a774.png

In [ ]:
wc -l *.count

b0eea28acf0f3a41b556a4adee2f443.png

看下每个文件的格式,查看前4行,第一列ensembl_gene_id,第二列read_count计数

In [ ]:
head -n 4 *.count

f0913a1d6235ee83ecf2bafb0b35b77.png

后4行

In [ ]:
tail -n 4 *.count

02ae9a798663033316fbef52abe38b5.png

8.1 合并表达矩阵并进行基因名注释(Rstudio中进行)

Tumor:SRR316214,SRR316215

Adjacent Normal Liver:SRR316212,SRR316213

In [5]:
options(stringsAsFactors = FALSE)
control1<-read.table("SRR3589959.count",sep = "\t",col.names = c("gene_id","control1"))
In [6]:
head(control1)
A data.frame: 6 × 2
gene_idcontrol1
<chr><int>
1ENSMUSG00000000001.4 1648
2ENSMUSG00000000003.15 0
3ENSMUSG00000000028.15 835
4ENSMUSG00000000031.16 65
5ENSMUSG00000000037.17 71
6ENSMUSG00000000049.11 0
In [8]:
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行是我们不需要的,可以删除

In [9]:
tail(control1)
A data.frame: 6 × 2
gene_idcontrol1
<chr><int>
56262ENSMUSG00000118659.1 1
56263__no_feature 12740305
56264__ambiguous 278925
56265__too_low_aQual 1099164
56266__not_aligned 498913
56267__alignment_not_unique13197674

进行整合

In [10]:
raw_count <- merge(merge(control1, control2, by="gene_id"), merge(treat1, treat2, by="gene_id"))
In [11]:
head(raw_count)
A data.frame: 6 × 5
gene_idcontrol1control2treat1treat2
<chr><int><int><int><int>
1__alignment_not_unique13197674164263932229671518776033
2__ambiguous 278925 344834 476611 374079
3__no_feature 12740305152224402249638618850217
4__not_aligned 498913 441735 684272 404241
5__too_low_aQual 1099164 1293547 1783744 1522713
6ENSMUSG00000000001.4 1648 2306 2941 2780
In [12]:
tail(raw_count)
A data.frame: 6 × 5
gene_idcontrol1control2treat1treat2
<chr><int><int><int><int>
56262ENSMUSG00000118654.1 0 0 0 0
56263ENSMUSG00000118655.111111710
56264ENSMUSG00000118656.1 0 0 0 0
56265ENSMUSG00000118657.1 0 0 0 0
56266ENSMUSG00000118658.1 0 0 0 0
56267ENSMUSG00000118659.1 1 4 6 1

删除前5行

In [14]:
raw_count_filt <- raw_count[-1:-5,]
head(raw_count_filt)
A data.frame: 6 × 5
gene_idcontrol1control2treat1treat2
<chr><int><int><int><int>
6ENSMUSG00000000001.4 1648230629412780
7ENSMUSG00000000003.15 0 0 0 0
8ENSMUSG00000000028.15 835 94913661051
9ENSMUSG00000000031.16 65 83 52 53
10ENSMUSG00000000037.17 71 53 94 66
11ENSMUSG00000000049.11 0 3 4 5

因为我们无法在EBI数据库上直接搜索找到ENSMUSG00000024045.5这样的基因,只能是ENSMUSG00000024045的整数,没有小数点,所以需要进一步替换为整数的形式。

In [16]:
# 第一步将匹配到的.以及后面的数字连续匹配并替换为空,并赋值给ENSEMBL
ENSEMBL <- gsub("\\.\\d*", "", raw_count_filt$gene_id) 
# 将ENSEMBL重新添加到raw_count_filt1矩阵
row.names(raw_count_filt) <- ENSEMBL
head(raw_count_filt)
A data.frame: 6 × 5
gene_idcontrol1control2treat1treat2
<chr><int><int><int><int>
ENSMUSG00000000001ENSMUSG00000000001.4 1648230629412780
ENSMUSG00000000003ENSMUSG00000000003.15 0 0 0 0
ENSMUSG00000000028ENSMUSG00000000028.15 835 94913661051
ENSMUSG00000000031ENSMUSG00000000031.16 65 83 52 53
ENSMUSG00000000037ENSMUSG00000000037.17 71 53 94 66
ENSMUSG00000000049ENSMUSG00000000049.11 0 3 4 5

保存count矩阵,方便后续进行下游分析

In [17]:
write.csv(raw_count_filt, file='raw_count_filt.csv')