RNAseq常规流程 |
您所在的位置：网站首页 › 测序比对率低会导致什么 › RNAseq常规流程 |
|
转录组分析是生信分析的一个基础和常见的分析流程，一般从下机数据开始，经过一系列的质控，基因组比对，差异分析等过程得到实验组与对照组（或者肿瘤中的转移组和原癌组）中差异表达的基因集，然后在基于该基因集结合我们的研究方向进行进一步的功能分析，信号通路分析，细胞间通讯，实验验证等得到一个可解释的实验结果。这篇文章主要就涉及到的基本流程和相应软件进行一个简单的介绍，后续更深入的知识点还需要小伙伴们自己去挖掘哦。 基本分析流程： 原始数据 -- Rawdata 质控 -- QC -- Trimmomatic，FastQC 基因组比对 -- Aligment -- STAR 比对结果量化 -- FeatureCount -- featurecount 两组间的基因表达差异分析 -- DEGs -- limma 差异基因的功能富集分析 -- GO,KEGG -- clusterProfiler 一. 测序数据下机后处理首先，我们拿到的原始序列文件就是fq.gz 文件 一般研究人员会把自己的样本（实验组/对照组）送到测序公司进行测序，然后测序公司会返回给我们一个fq.gz文件，如下. 每个样本都有两条链的测序数据 *_1.fq.gz 和 *_2.fq.gz ， 我们要把他们单独放在一个文件夹下作为一个样本进行后续分析。 fq.gz.png 1.下载fastq.gz rawdatawget url 2.每个样本创建一个以样本名命名的文件夹（eg RNA-1）, 并将相应的._1.fq.gz/._2.fq.gz 放入文件夹（eg RNA-1），用作后续QC处理。 mkdir sample1 sample2 sample3 control1 control2 control3 for samp in {sample1 sample2 sample3 control1 control2 control3} do mv "../rawdata/*"$samp"*_1.fq.gz" $samp mv "../rawdata/*"$samp"*_2.fq.gz" $samp done 二. 低质量数据过滤和质量评估QC -- Trimmomatic/FastQC 1. TrimmomaticTrimmomatic: A flexible read trimming tool for Illumina NGS data 因为公司给的数据里面，是没有去掉接头的，所以需要跑之前先把接头去掉，过滤质量后再用。 使用Trimmomatic(或者Trim Galore)切除数据中的接头序列和低质量序列，专门清洗illumina测序数据的常用工具 优点： 1、可直接对.fq.gz压缩文件操作 2、结果默认生成在同文件的目录下，不生成文件夹，直接生成四个*.fq.gz文件。若需要有文件夹，需要额外建立文件夹的指令。参数-o,为文件夹（文件） trim_result.png Trimmomatic的执行文件是一个Java文件 ### 切除接头序列 java -jar /home/Software/Trimmomatic-0.39/trimmomatic-0.39.jar \ PE \ # rawdata为双端测序 -phred33 \ # Fastq文件的质量值格式 -threads 6 \ # 6 线程 -trimlog Trimmomatic.log \ # 日志文件保存位置 $read1 $read2 \ # *_1.fq.gz, *_2.fq.gz 需要过滤的Fastq文件 $trim_read1 \ # *_1.clean.fq.gz 过滤后的Fastq文件 output_unpaired_1.fq.gz \ # read1的5‘->3'和read2的3‘->5'未能匹配上的Fastq文件 $trim_read2 \ output_unpaired_2.fq.gz \ ILLUMINACLIP:/home/Software/Trimmomatic-0.39/adapters/TruSeq3-PE-2.fa:2:30:10:1:TRUE \ # 去除illumina测序平台下的TruSeq3接头序列，不知道可以问测序公司。 # 2:30:10即表示，在比对接头序列时允许有两个位置的碱基发生错配，双端测序的两条reads与接头序列匹配率超过30%的话，就会被切除掉，单条reads如果与接头序列的匹配率超过10%，也会被切除掉。 LEADING:0 \ # 切除reads 5’端质量值低于0的碱基, 即不对5’端剪切 TRAILING:25 \ # 切除reads 3’端质量值低于25的碱基 SLIDINGWINDOW:50:25 \ #以50bp为窗口进行滑窗统计，切除碱基平均质量低于25的窗口及之后的序列。 MINLEN:50 \ #去除过滤后长度低于50的reads HEADCROP:0 \ # 切除reads开头的0个碱基 2>trim.log 1>trim.err # 0 表示stdin标准输入；1 表示stdout标准输出；2 表示stderr标准错误使用参数 PE/SE rawdata为 Paired-End，还是 Single-End的reads，其输入和输出参数稍有不同。 -threads 设置多线程运行数 -phred33 设置碱基的质量格式，可选pred64 ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 切除adapter序列。参数后面分别接adapter序列的fasta文件：允许的最大mismatch数：palindrome模式下匹配碱基数阈值：simple模式下的匹配碱基数阈值。 LEADING:3 切除首端碱基质量小于3的碱基 TRAILING:3 切除尾端碱基质量小于3的碱基 SLIDINGWINDOW:4:15 Perform a sliding window trimming。Windows的size是4个碱基，其平均碱基质量小于15，则切除。 MINLEN:50 最小的reads长度 CROP: 保留reads到指定的长度 HEADCROP: 在reads的首端切除指定的长度 TOPHRED33 将碱基质量转换为pred33格式 TOPHRED64 将碱基质量转换为pred64格式 2.FastQC (测序数据质控) FastQC 是高通量测序数据的高级质控工具，输入FastQ，SAM，BAM文件，输出对测序数据评估的网页报告。 一般在Trimmomatic之后，用FastQC软件对预处理数据进行质量控制分析。 优点： 1、可直接对.fq.gz压缩文件操作 2、结果默认路径生成在同文件的目录下，生成的文件夹为XX_fastqc,里面有fastqc_data.txt , fastqc_report.html等重要文件 fastqc $trim_read1 $trim_read2 fqc.trim.png fqc_detail.trim.png 3、可生成在指定的文件夹（应该说是路径吧），需要额外建立文件夹，用-o参数，为文件夹（文件） fastqc -t 12 -o out_path sample1_1.fq sample1_2.fq 使用说明 -t --threads:线程数 -o --outdir:输出路径 --extract：结果文件解压缩 --noextract：结果文件压缩 -f --format:输入文件格式.支持bam,sam,fastq文件格式结果解读 我们对Trim过滤之后的序列进行质量评估，评估结果如下： FastQC 会生成一个fastqc_report.html的结果报告，下面是软件对质控结果进行判断：绿色代表PASS；黄色代表WARN；红色代表FAIL image.png简单说一下报告中常见的几张图吧 每个碱基位置的质量分数分布图 image.png横轴为read长度，纵轴为质量得分。说明trim后的Reads质量较高，可以用于后续分析。 一般要求所有位置的10%分位数大于20，即大于最多允许该位置10%的序列低于Q20。绿色区域说明数据质量较高。一般但随着测序长度的加深，碱基的错赔率会增加。 平均质量分数分布图 image.png横轴表示Q值，纵轴表示每个值对应的read数目，当测序结果主要集中在高分中，证明测序质量良好 每个碱基位置的碱基N百分比 image.pngN占比越低，数据质量越高 序列的核酸组成偏好性 image.png正常情况下，我们期望的是A=T=C=G，GC含量较高时，经常会出现A=T,C=G的情况。 每条reads的GC占比 image.pngGC含量-GC含量对应的read数。蓝色为经验分布给出的理论值，红色是测序数据的真实值，但两条线接近重合时，说明数据质量较高；当红色出现双峰可能是样品被污染或不纯。 重复reads统计 image.png横轴表示重复的次数，纵轴表示重复的reads的数目 其他更多信息可以去官网自行查看。 三. Clean reads的基因组比对 Aligment -- STAR针对每个样本，利用 STAR软件 将预处理序列与测序物种的参考基因组序列进行序列比对。 1. 比对原理将QC后的reads比对到参考基因组上，一般分为有参考基因组比对和无参考基因组比对。我们常用的就是有参考基因组的比对。 aligment.png 如果数据库中存在该物种的物种基因组（fasta）及其基因组注释文件（gff, gtf），那么，我们就可以用该基因组信息作为参考基因组进行clean reads序列比对，便可以知道基因组的哪些位置转录了RNA序列，通过对reads进行计数，而得到基因表达的定量数据。 如果没有可用的参考基因组及其注释信息，则需要考虑无参考基因组的比对方法进行比对。首先先将reads 进行组装拼接，得到从头组装的组装转录本，然后再将clean reads比对到组装转录本上得到比对结果。在参考基因组比对的时候，会优先考虑大片段比对，如果比对不上，再考虑将片段切割成不同区段进行比对。 RNAseq 常用的比对工具有：Hisat2, STAR 等 2. STAR 比对 STAR --runThreadN 20 \ # 线程数 --genomeDir /home/database/ref \ # 参考基因组 --readFilesCommand gunzip -c \ # *.gz 文件解压缩 --readFilesIn *_1.fq.gz *_2.fq.gz \ #输入文件，空格隔开 --outSAMtype BAM SortedByCoordinate \ # 输出bam文件（默认输出sam文件） --outFileNamePrefix Aligment \ # 输出文件的位置和前缀结果解读 比对结果：Aligned.out.sam或者Aligned.out.bam 比对完以后的统计信息： Log.final.out 剪切的信息：SJ.out.tabBAM/SAM文件 bam是sam的压缩格式，在内容上是一样的，linux下可使用samtools, bamtools 进行查看与转换。 bam文件分为2个部分，header和record. header,以@开头，包含： @HD 版本信息； @SQ 比对参考序列信息； @RG 测序平台信息； @PG 序列比对使用的软件信息record为基因组比对信息，每一行信息为： reads ID 比对到基因组上的信息 比对到的染色体 比对上的位置信息 比对质量值 碱基组成 碱基质量值等 bam.png比对率统计 mapcount.png总比对率和唯一比对率是我们关注的，尤其是总比对率。 总比对率低，一般可分为2种情况，总比对率低于10% ，可能是弄错了物种来源，比如鼠源细胞当成了人源的；总比对率低于60%，说明有外源基因组的污染，最常见的就是细胞培养的时候受到支原体污染。比对区域分布统计 image.png比对reads在每条染色体上的密度分布情况 image.png可以通过IGV软件可视化查看BAM文件信息 具体以后再详细展开…… 四. 比对结果的外显子（基因）表达计数 FeatureCount -- featurecount在高通量测序分析中用于下游分析的关键信息是比对到每个genomic feature（外显子的表达量、featurecounts是一款使用于RNA-seq和DNA-seq的read summarization工具，应用了高效率的染色体哈希算法和feature区块技术因等）中的read数目。 featurecounts 速度快，需要的内存空间少，同时可以用于单端和双端的数据 通过STAR等软件比对得到的BAM文件，利用计数软件（如 featurecount, HTseq, stringtie），获取各基因的 reads count。 假设前提：认为比对到基因的reads 数与该基因的表达量成正相关。 readscount.pngfeaturecount的调用方法 /home/Software/featureCounts \ -T 4 \ # 线程数目，1~32 -a $gtf \ # 参考gtf文件 -o read.count \ # 输出文件的名字，输出文件的内容为read 的统计数目 -p \ # 只能用在paired-end的情况中，会统计fragment而不统计read (-p -B -C 同时使用、不使用) -B \ # 在-p选择的条件下，只有两端read都比对上的fragment才会被统计 -C \ # 如果-C被设置，那融合的fragment就不会被计数，这个只有在-p被设置的条件下使用 -t exon \ # 设置feature-type，-t指定的必须是gtf中有的feature，同时read只有落到这些feature上才会被统计到，默认是“exon” -g gene_name \ # 默认为gene_id，注意！选择gtf中提供的id identifier！！！ -s 2 \ # 2条链比对 Align.sorted.bam # 输入的bam/sam文件，支持多个文件输入得到read.count数据之后，还需要进行进一步处理才能够进行分析。 由于实验中每个样本的上样量不同，或者基因的长度,测序深度不同等导致基因count本身就存在着较大差异，所以，我们是不能够直接只用count 数据进行后续的基因表达差异分析的。因此我们需要对 count 数据进行标准化。 countToFpkm |
今日新闻 |
推荐新闻 |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |