RNA |
您所在的位置：网站首页 › 基因富集分析软件 › RNA |
RNA
|
今天分享的学习笔记是一套转录组分析简单流程，适用于初学者入门阅读，从原始测序数据开始，经过质控、序列比对、定量表达、差异表达、功能富集等一系列分析步骤，最终获得基因表达信息。本文所有数据都经过特殊修改，仅供学习参考使用。 转录组是在特定时空条件下细胞中基因转录表达产物，广义的转录组包括信使RNA，核糖体RNA，转运RNA及非编码RNA，狭义上是指所有mRNA的集合，转录组分析能够获得不同基因的表达情况。 1. 数据来源假设有两个不同组织（PR和SR），每个组织各区三个样本，一共六个样本，利用illumina平台进行转录组测序，得到双端测序数据。数据原始格式为.fq,共有12条测序数据文件（每个样本产生两条） PR1_2.fq PR2_2.fq PR3_2.fq SR1_2.fq SR2_2.fq SR3_2.fq PR2_1.fq PR3_1.fq SR1_1.fq SR2_1.fq SR3_1.fq创建工作文件夹，并新建子步骤目录：00_trainingRawdata 01_trimmomaticFiltering 02_hisat2Mapping 03_featurecountsQuatification 04_DESeq2DEGanalysis mkdir RNA-Seq_Practice cd RNA-Seq_Practice mkdir 00_trainingRawdata 01_trimmomaticFiltering 02_hisat2Mapping 03_featurecountsQuatification 04_DESeq2DEGanalysis #批量创建工作文件夹 2. 测序数据质量评估利用fastQC软件对获得的fastq序列文件进行质量分析，生成html格式的结果报告，其中含有各项指标，以下用PR1样本为例。 fastqc *.fq #批量对fq后缀文件运行fastqc程序 输出结果：PR1_1_fastqc.html Filename PR1_1.fq File type Conventional base calls Encoding Illumina 1.5 Total Sequences 105300 #总序列数 Sequences flagged as poor quality 0 Sequence length 90 #序列长度 %GC 52 #GC碱基含量质量评估报告，使用浏览器打开：
利用Trimmomatic软件除去序列文件中的接头（adapter），并对碱基进行合适的修改，然后对碱基进行修剪，对低质量的序列进行过滤。 time java -jar trimmomatic-0.39.jar PE #trimmomatic是依赖java环境运行 -threads 1 -phred33 PR1_1.fq PR1_2.fq -summary ../01_trimmomaticFiltering/PR1.summary -baseout ../01_trimmomaticFiltering/PR1 LEADING:3 TRAILING:3 SLIDINGWINDOW:5:20 HEADCROP:13 MINLEN:36运行程序对序列文件中低质量的序列进行过滤，将输出结果存储到01_trimmomaticFiltering，每一个样本序列会生成如下输出文件，summary文件包含过滤的总结信息。 -rw-r--r-- 1 19M Nov 28 14:29 PR1_1P -rw-r--r-- 1 0 Nov 28 14:29 PR1_1U -rw-r--r-- 1 19M Nov 28 14:29 PR1_2P -rw-r--r-- 1 0 Nov 28 14:29 PR1_2U -rw-r--r-- 1 282 Nov 28 14:29 PR1.summary打开其中一个PR1.summary文件进行查看，其中Input Read Pairs表示过滤之前的reads数，Both Surviving Reads表示过滤之后的reads数。 $ cat PR1.summary #查看总结文件 Input Read Pairs: 102300 Both Surviving Reads: 102300 Both Surviving Read Percent: 100.00 Forward Only Surviving Reads: 0 Forward Only Surviving Read Percent: 0.00 Reverse Only Surviving Reads: 0 Reverse Only Surviving Read Percent: 0.00 Dropped Reads: 0 Dropped Read Percent: 0.00 4. 比对到参考基因组利用hisat2软件，将fasta序列比对到参考基因组。首先需要构建索引文件，下载或者拷贝参考基因组序列文件Ref和基因组注释文件，通过hisat2-build命令构建索引文件。 cd xx/RNA-Seq_Practice cp -r xx/Ref . cd Ref gunzip -c chr.fa.gz > chr.fa #解压参考基因组 time hisat2-build -p 1 chr.fa Chr #建立索引文件完成上述步骤后，将过滤得到的reads比对到参考基因组上。输入文件为两个fasta序列文件，将运行过程中的输出提示重定向到log文件。 -S设置输出文件为sam-x设置参考基因组-p设置运行的线程数量 time hisat2 -p 1 #线程数为1 -x Ref/Chr #参考基因组文件路径 -1 01_trimmomaticFiltering/PR1_1P #输入文件 -2 01_trimmomaticFiltering/PR1_2P #输入文件 -S 02_hisat2Mapping/PR1.sam #输出文件sam格式 --new-summary 1>02_hisat2Mapping/PR1_hisat2Mapping.log 2>&1 #输出日志得到的日志文件中包含比对成功的reads数量和比对率等信息： HISAT2 summary stats: Total pairs: 102300 Aligned concordantly or discordantly 0 time: 94209 (92.09%) Aligned concordantly 1 time: 7613 (7.44%) Aligned concordantly >1 times: 293 (0.29%) Aligned discordantly 1 time: 185 (0.18%) Total unpaired reads: 188418 Aligned 0 time: 186684 (99.08%) Aligned 1 time: 1575 (0.84%) Aligned >1 times: 159 (0.08%) Overall alignment rate: 8.76% #总比对率8.76%上述步骤完成后，会在当前目录下生成多个sam格式文件，SAM（sequence Alignment/mapping)格式是高通量测序中存放比对数据的标准格式。另外，bam是sam的二进制格式，可以减小sam文件的大小，因此利用samtools对sam进一步处理，得到bam文件，以下用PR1为例，其他样本按相同方式处理。 time samtools view -bS 02_hisat2Mapping/PR1.sam > 02_hisat2Mapping/PR1.bam time samtools sort 02_hisat2Mapping/PR1.bam > 02_hisat2Mapping/PR1.srt.bam 5. 基因表达定量分析利用featuresCounts软件，对reads进行精确的计数，最后将所有样本的reads数合并为一个文件，将RNA-Seq_Practice_countstable文件导出，根据FPKM和TPM的计算公式定量 每个基因在每个样本中的表达量。利用网站（https://www.ncbi.nlm.nih.gov/gene）将基因ID转化为Gene description，从而了解其功能和相关信息。 time featureCounts -a RefMaize_AGPv4/Zea_mays.AGPv4.32.gtf.gz #基因组注释文件路径 -T 1 -p --countReadPairs -g gene_id -t exon -o 03featurecountsQuatification/PR1 02hisat2Mapping/PR1.srt.bam以PR1为例，将多余的信息除去，只保留基因ID和reads数，得到count文件，然后同样步骤处理其他样本，最后将所有的reads定量数据合并为一个countstable文件。 cat PR1 | cut -f 1,7 > PR1.count paste PR1.count PR2.count PR3.count SR1.count SR2.count SR3.count > countstable #合并不同样本的定量信息到一个文件 awk '{$3="";$5="";$7="";$9="";$11="";print $0}' countstable > RNA-Seq_Practice_countstable #提取指定列到新文件将生成的RNA-Seq_Practice_countstable保存到本地，然后计算FPKM和TPM值，在R语言中进行相关计算。 FPKM（Fragments Per Kilobase of exon model per Million mapped fragments）表示每千个碱基的转录每百万映射读取的fragments，该方法是利用每个样本的总fragments数进行校正。优点是消除样本间和基因间的差异，但是该方法主要是计算的结果可能与真实表达水平有偏差。 
TPM (transcript per million) 表示每百万转录本中来自于某基因的转录本数目，该方法先对每个基因的reads数用基因的长度进行校正，然后再利用校正后的reads数与校正后该样本所有的reads数求商。优点是消除了外显子长度造成的差异和样本间测序总reads counts不同造成的差异，缺点为该法不是采用比对到基因组上的总reads counts，所以有时候不够精确。 
具体计算方法如下所述： > df rownames(df) df colnames(df) dim(df); names(df) #查看df的数据结构和变量名称是否正确 [1] 41768 6 [1] "PR1" "PR2" "PR3" "SR1" "SR2" "SR3" > head(df[,1:6],2) #输出前两行数据进行预览 PR1 PR2 PR3 SR1 SR2 #PR和SR共六列数据 Zm00001d027 0 0 0 0 0 Zm00001d021 0 0 0 0 0 featureCounts_meta 8. KEGG分析在R语言（4.2.2）中使用clusterProfiler、KEGG.db、ggplot2三个R包进行KEGG富集分析，先构建一个OrgDb本地索引库，然后利用enrichKEGG函数进行富集分析，利用在线网站(https://www.maizegdb.org/gene_center/gene)将基因ID转换为GenBank Gene类型，保存为文件gene.txt，输出文件为pathway数据表gene.pathway.csv 和KEGG富集结果gene.pathway .pdf remotes::install_github("YuLab-SMU/createKEGGdb",force = TRUE) createKEGGdb::create_kegg_db('zma') #首先使用该软件制作玉米KEGG富集包 library(clusterProfiler) library(ggplot2) # 绘图需用R包 library(KEGG.db) # 该步骤需要手动在Rstudio中安装KEGG.db包 file="gene.txt" #输入文件，包含18条差异表达的基因信息 gene head(gene,2) GenBank.Gene Zm00001d0218 1036647 Zm00001d0234 1002408 kk 9. 差异表达基因火山图通过ggplot2软件绘制基因差异表达火山图，使用DESeq2软件生成的matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results为输入文件，该文件内geneid、sampleA、sampleB、baseMeanA、baseMeanB、baseMean、log2FoldChange、lfcSE、stat、pvalue、padj等信息，将其导入R语言中，提取1,6,7,8,9,10,11列，然后设置数据的行列名称 library(ggplot2) # 载入R包 res |
【本文地址】
今日新闻 |
推荐新闻 |