RNAseq analysis¶
上游分析¶
RNA数据的获取¶我们通常分析的是DNA转录后所对应的RNA链,一般从NCBI数据库获取SRA数据来进行RNAseq分析.因此在获取SRA编号后可在服务器中下载.下载方式有以下几种：
aspera 工具下载
wget, curl 命令直接下载
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]
--quality
--phred33 #设定记分方式，代表Q+33=ASCII码的方式来记分方式。这是默认值。
--paired # 对于双端结果，一对reads中若一个read因为质量或其他原因被抛弃，则对应的另一个read也抛弃.
--output_dir #输出目录，需确保路径存在并可以访问
--length #设定长度阈值，小于此长度会被抛弃.这里测序长度是100我设定来75，感觉有点浪费
--strency #设定可以忍受的前后adapter重叠的碱基数，默认是1.不是很明白这个参数的意义
-e
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下载.
GFF(general feature format)：用于基因组注释.
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
参数说明：
-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
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 [ ]: