RNAseq analysis

上游分析

RNA数据的获取

我们通常分析的是DNA转录后所对应的RNA链,一般从NCBI数据库获取SRA数据来进行RNAseq分析.因此在获取SRA编号后可在服务器中下载.下载方式有以下几种：

  1. aspera 工具下载
  2. wget, curl 命令直接下载
  3. NCBI官方的 SRA Toolkit 进行下载 https://zhuanlan.zhihu.com/p/89024212

SRA数据转换成fastq文件

获取sra数据后在服务器中用fastq-dump 命令对其进行转换，转换后得到fastq文件,该文件可进行RNAseq的上游分析.

fastq-dump --gzip --split-3 -O path -A SRR1039508

md5验证转换的fastq.gz文件是否完整

md5sum *gz md5.txt #给自己的文件生成md5值 md5sum -c md5.txt #比对已有的md5值 若结果均显示OK,则表示该数据文件完整

fastqc检测测序文件质量

fastqc -o outputdir SRR1039508_1.fastq.gz SRR1039508_2.fastq.gz

查看QC文件

质控之后生成zip文件以及html文件,可查看网页报告，若RNA为双链，会生成两个相应的文件.https://zhuanlan.zhihu.com/p/88655260

multiqc 质量报告（合并多个fastq报告）

multiqc path/*zip

rawdata的过滤和清除不可信数据---trim_galore

trim_galore：可以处理illumina，nextera3，smallRNA测序平台的双端和单端数据，包括去除adapter和低质量reads.

In [ ]:
trim_galore [options] <filename>
--quality<int>  #设定phred quality阈值。默认20（99%的read质量），如果测序深度较深，可以设定25
--phred33       #设定记分方式，代表Q+33=ASCII码的方式来记分方式。这是默认值。
--paired          # 对于双端结果，一对reads中若一个read因为质量或其他原因被抛弃，则对应的另一个read也抛弃.
--output_dir   #输出目录，需确保路径存在并可以访问
--length        #设定长度阈值，小于此长度会被抛弃.这里测序长度是100我设定来75，感觉有点浪费
--strency     #设定可以忍受的前后adapter重叠的碱基数，默认是1.不是很明白这个参数的意义
-e<ERROR rate>  #设定默认质量控制数，默认是0.1，即ERROR rate大于10%的read 会被舍弃，如果添加来--paired参数则会舍弃一对reads
<filename>  #如果是采用illumina双端测序的测序文件，应该同时输入两个文件.

RNA单链

~/TrimGalore-0.4.5/trim_galore -o path -U SRRxxxxxxx

RNA双链

~/TrimGalore-0.4.5/trim_galore -o path --paried SRR1039508_1.fastq.gz SRR1039508_2.fastq.gz

In [ ]:
nohup ~/TrimGalore-0.4.5/trim_galore -o path --paried SRR1039508_1.fastq.gz SRR1039508_2.fastq.gz>ssr.out&(挂起)

整理后数据的质控

对过滤后对文件进行质量分析。观察过滤结果。同样使用fastqc和multiqc两个软件进行质量分析.过滤后数据写入 SRR1039508_1_val_1.fq.gz&SRR1039508_2_val_2.fq.gz中.

In [ ]:
fastqc path/*gz -o path

使用hisat2比对回帖

1.建立索引

如果自己建立索引文件的话，需要消耗大量时间，通常人类基因组和一些动植物的基因组的索引文件可以在NCBI、Ensembl、UCSC、GeneCode下载

2.基因注释文件GTF/GFF的下载

基因注释文件也可在NCBI、Ensembl、UCSC、GeneCode下载.

  1. GFF(general feature format)：用于基因组注释.
  2. GTF(gene transfer format)：用于对基因的注释. https://blog.csdn.net/shandg_lxy/article/details/89182341 https://zhuanlan.zhihu.com/p/79631226
In [ ]:
#下载参考基因组
wget https://ftp.ncbi.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.28_GRCh38.p13/GCA_000001405.28_GRCh38.p13_genomic.fna.gz 
gzip -d GCA_000001405.28_GRCh38.p13_genomic.fna.gz 
 
#下载gff注释文件
wget https://ftp.ncbi.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.28_GRCh38.p13/GCA_000001405.28_GRCh38.p13_genomic.gff.gz
gzip -d GCA_000001405.28_GRCh38.p13_genomic.gff.gz

hisat2回帖

In [ ]:
hisat2  -p 6    -x <dir of index of genome>  -1  SRR1039508_1_val_1.fq.gz   -2  SRR1039508_2_val_2.fq.gz  -S  tem.hisat2.sam
参数说明：

-p #多线程数 -x #参考基因组索引文件目录和前缀 -1 #双端测序中一端测序文件 -2 #同上 -S #输出的sam文件

说明：在比对过程中，hisat会自动将双端测序匹配同一reads并在基因组中比对，最后两个双端测序生成一个sam文件。比对回帖过程需要消耗大量时间和电脑运行速度和硬盘存储空间。5G左右fastq文件比对回帖过程消耗大概一个小时，生成了17G的sam格式文件。回帖完成会生成一个回帖报告。

samtools软件进行格式转换

SAM文件和BAM文件 samtools 是针对比对回帖的结果——sam和bam格式文件的进一步分析使用的软件。sam格式文件由于体量过大，一般都是使用bam文件来进行存储。由于bam文件是二进制存储所以文件大小比sam格式文件小许多，大约是sam格式体积的1/6 。 https://www.jianshu.com/p/d978fddb9a45

In [ ]:
samtools view -bS seq.sam  > seq.bam  #文件格式转换
samtools sort seq.bam  seq_sorted.bam  ##将bam文件排序
samtools index seq_sorted.bam  #对排序后对bam文件索引生成bai格式文件，用于快速随机处理.
samtools view -h s.bam|less -S
samtools view s.bam|less -S
# 提取chr1染色体，生成只有chr1的bam文件
samtools view -h -b s.bam chr1 >s.chr1.bam
samtools view -bt ref_list.txt -o aln.bam aln.sam.gz
In [ ]:
#本文例子的linux操作
samtools view -bS tem.hisat2.sam >tem.bam
samtools sort tem.bam tem.sorted.bam
samtools index tem.sorted.bam tem.sorted.bam

到这一步一个回帖到基因组对RNA-seq文件构建完成.

对回帖bam文件进行质量评估

samtools falgstat ：统计bam文件中比对flag信息，然后输出比对结果.

In [ ]:
samtools flagstat tem.sorted.bam > SRR1039508.sorted.flagstat
gtf<-list.files("./",pattern=".gtf$",full.names=T)
fetureCounts(bam.files,annot.ext=gtf,isGTFAnnotationFile=T)

计数（Count）

计算RNA-seq测序reads对在基因组中对比对深度. 计数工具：subread中的featureCounts.

In [ ]:
featureCounts -T 6 -t exon -g gene_id -a <gencode.gtf> -o fc.txt tem.sorted.bam#count
cut -f 1,7 fc.txt |grep -v '^#'>SRR1039508.counts.txt#单独生成counts文件
cat SRR1039508.counts.txt #查看txt

-g # 注释文件中提取对Meta-feature 默认是gene_id -t # 提取注释文件中的Meta-feature 默认是 exon -p #参数是针对paired-end 数据 -a #输入GTF/GFF 注释文件 -o #输出文件

这一步获得的fc.txt可直接用于下游分析,即可在R的环境下构建表达矩阵后，进行差异性分析、富集分析、聚类分析.

下游分析

In [ ]: