RASflow原理理解、项目部署及运行

Quasimodo 生信菜鸟团 2023-06-07 12:30 发表于湖北

RASflow原理理解、项目部署及运行--for GSE103584


本周转录组专辑将转载一位老师使用RASflow一体化、流程化分析的笔记，笔记内容非常详实，让我这种散装户受益匪浅

（本文所有引用部分为我个人的一些记录，读者需要注意区分于正文，感兴趣的小伙伴可以直接点击“阅读原文”）

RASflow流程在生信技能树的公众号之前的推文中也提到过：分享你的NGS数据分析流程也能发文章哦

感兴趣的小伙伴可以看看这篇2020年的推文和原始文献

流程整体的理解

流程重点

一、两次python main.py，QC yes做质控 （根据QC决定是否行数据修剪trim QC） no 质控后流程

  • 若QC合格，则文中推荐trim=no，因为数据修剪可能会丢失信息
  • If the quality of the reads is good enough, it is recommended that trimming should not be performed since it would lead to loss of information; but if the quality is low,trimming is suggested to improve the quality
  • trimming is performed using the tool Trim Galore https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/

二、map to genome or transciptome?

  • transcriptome：则pseudo alignment and quantify with salmon【2017nature methods】 https://www.nature.com/articles/nmeth.4197
  • genome：alignment with hisat2 and count with featureCounts/htseq-count

三、注释文件的选择与下载：不同的注释文件可以影响比对准确性

（一）注释文件包括：①基因组序列【fasta 】②基因组注释文件【gtt/gtf】
  • ①基因组序列【fasta 】就是那3×109的TCGA碱基序列
  • ②基因组注释文件【gtt/gtf】就是对碱基序列区段的注释：哪个基因，发挥什么功能
（二）不同机构对基因组注释文件命名有差异

一文读懂参考基因组和基因组注释+最全下载方法_白墨石的博客-CSDN博客 https://blog.csdn.net/u011262253/article/details/117486244

图片
（三）注释文件选择大原则
  • 强调基因表达估计的可重复性和稳健：优选较为简洁的、NCBI中的RefSeq【提供了便于识别和遴选的下载链接】的fasta和gtf

https://www.ncbi.nlm.nih.gov/datasets/taxonomy/9606/

2022-Impact of gene annotation choice on the RNAseq.pdf(1.7 MB)

  • 强调进行更具探索性，需要更全面的注释：Ensembl http://www.ensembl.org/info/data/ftp/index.html【如何选择，见下】
  • GeneCode https://www.gencodegenes.org/ 有对参考序列和注释文件的描述，适合阅读和研判
（四）NCBI中的RefSeq,GenBank及all包含的基因序列及注释文件的来源

https://support.nlm.nih.gov/knowledgebase/article/KA-03339/en-us

要访问当前和正在更新的基因组组装数据，请使用NCBI基因组FTP网站上的以下三个目录：genbank，refseq和all。

  • genbank是原始基因组组装数据的目录，包含测序中心或个人研究者提交给GenBank或国际核苷酸序列数据库协作组织（INSDC）的组装基因组序列和相关注释（如果有）。如果您想获取所有提交的基因组组装，并且主要关注的不是访问基因组注释，请使用此目录。该目录按分类群组织，您可以直接浏览它。
  • refseq是NCBI衍生的基因组组装数据目录，包含NCBI RefSeq工作人员从主要INSDC数据中选择的已组装基因组。如果您对高质量且经常维护的注释数据感兴趣，请使用refseq目录。RefSeq基因组组装的序列是与相应INSDC组装中存在的序列相同的副本。在某些情况下，复制可能不完全相同，因为RefSeq工作人员可能会（1）删除序列的较小片段（称为contigs）或报告的污染物，或者（2）将非核基因组序列（例如线粒体）添加到组装中。要查找用于创建RefSeq组装的主要GenBank（INSDC）组装，请使用装配报告文件。所有RefSeq基因组组装都具有注释，RefSeq工作人员通过NCBI原核或真核基因组注释管道提供或从主要记录传播。refseq目录中存在的基因组组装数量小于genbank目录中的数量。该目录按分类群组织，您可以直接浏览它。
  • all是将genbank和refseq目录内容结合在一起的目录。它由两个主根目录组成：GCA用于GenBank和GCF用于RefSeq组装数据。每个根目录都分成一个子目录层次结构，该层次结构遵循配件访问编号中数字的模式。例如，从GCF根开始，“000”目录将包含配件访问编号前三位数字为0、0、0的那些配件的子目录。“000/001/405”路径以一系列GCF_000001405目录结束，这些目录的名称反映了各个配件版本和名称（GCF_000001405.36_GRCh38.p10用于名为GRCh38.p10的人类参考配件版本36）。
  • NCBI基因组FTP网站上的所有其他目录都是旧目录，我们将逐步将它们存档。如果您正在使用这些目录之一，请注意它们的更新日期以确保您获取当前数据。如果发现缺少目录，请检查它是否已移动到存档目录中，您也会在Genomes FTP网站上找到它。阅读有关FTP基因组站点结构的更多信息，并了解有关网站重新组织、内容、文件格式、下载说明和未来计划的详细信息。
（五）Ensembl各个FASTA注释文件的含义

FASTA是一种常见的生物信息学文件格式，用于表示核酸序列或蛋白质序列。每个FASTA文件都由一系列的条目组成，每个条目以一个以">"字符开始的描述行，后面是序列行。

  1. DNA (FASTA)：这个文件包含了参考基因组的完整DNA序列。它包括了所有的染色体，以及任何已知的线粒体和质粒序列。
  2. cDNA (FASTA)：cDNA是指互补DNA，它是通过反转录过程从mRNA模板生成的。cDNA FASTA文件包含了所有已知基因的cDNA序列，它可以用于研究基因的表达。
  3. CDS (FASTA)：CDS是指编码序列（Coding Sequence），它是基因的那部分序列，可以被翻译成蛋白质。CDS FASTA文件包含了所有已知基因的编码序列。
  4. ncRNA (FASTA)：ncRNA是非编码RNA的缩写，它是指那些不编码蛋白质，但在其他生物过程中发挥重要作用的RNA。ncRNA FASTA文件包含了所有已知的非编码RNA序列。
  5. Protein sequence (FASTA)：这个文件包含了所有已知基因的蛋白质序列。它可以用于研究蛋白质的结构和功能。
（六）salmon通过pseudo-align进行transcript quantify:Ensembl 选 cDNA FASTA，此外，无需基因注释文件【我临时忘了，没想到chatgpt居然知道】

我们在往期推文MetaSRA：SRA的标准化元数据！中提到

refine.bio是一个存储库，存储来自公开来源的经过统一处理和规范化的转录组数据，它有自己处理不同来源SRA数据获取表达矩阵的pipelines

使用Salmon和tximport来处理refine.bio中的所有RNA-seq数据

Salmon是一种无比对的方法，用于从RNA-seq数据中估计转录物丰度。我们在准映射模式中使用它，这比基于比对的方法快得多，需要我们建立Salmon转录组索引。

Salmon的定量是在转录水平上。为了更好地与refine.bio中包含的微阵列数据整合，我们用tximport将转录物水平的信息总结到基因水平。

我们的tximport实现生成“lengthScaledTPM”，这是通过使用样本之间的平均转录长度和文库大小缩放TPM而生成的基因水平计数量表值。请注意，tximport是在实验级别应用的，而不是应用于单个样本。

如果你正在使用Salmon进行转录本定量，你应该选择cDNA (FASTA)文件。Salmon工作在转录本水平，需要cDNA序列来构建它的索引。这个cDNA FASTA文件包含了所有已知基因的cDNA序列，它可以用于研究基因的表达。

  • Salmon使用所谓的"pseudo-alignment"技术，这是一种可以快速精确地确定RNA-seq reads可能来自转录本哪一部分的方法。然后，Salmon利用这些信息计算每个转录本的表达量
图片
  • 在你下载了cDNA FASTA文件后，你可以使用Salmon的index命令来创建索引，然后使用quant命令来定量转录本。具体的命令可能像这样：
salmon index -t transcripts.fasta -i transcripts_index
salmon quant -i transcripts_index -l A -1 reads1.fq -2 reads2.fq -o quant

#在这里，transcripts.fasta是你下载的cDNA FASTA文件，
#transcripts_index是你想要创建的索引的名称，reads1.fq和reads2.fq是你的RNA-seq reads
#（如果你的数据是双端的话），而quant是你想要存储定量结果的目录。
# 注意，你需要根据你的实际情况来调整这些命令。例如，如果你的数据是单端的，
# 你只需要提供一个reads文件。此外，-l A选项告诉Salmon自动推断reads的库类型，
# 但如果你知道你的库类型，你可以直接指定它。

四、DEA：Deseq2/edgeR

（一）map to trancriptome：transcript-level DEA(Deseq2) ± gene-level DEA
  • transcript-level DEA：更推荐Deseq2包

    • DESeq2 is recommended for transcriptome-based and DEA because it provides specially a function DESeqDataSetFromTximport to import output from package tximport[3] which is used in transcriptome based pipeline as explained--《RASflow--Tutorial》
  • gene-level DEA：set GENE LEVEL为TRUE并依照tutorial作相应配置
（二）map to genome：DEA(Deseq2/edgeR均可)
（三）Deseq2 VS edgeR
图片

流程图

图片

正式开始

项目部署及运行

一、在服务器上通过docker进行环境部署

  1. docker安装，登陆并赋予sudo权力 https://www.yuque.com/duxiuju/upg5z6/by359d5wt3xl472k
  2. 通过docker拉取RASflow：docker container run --name=rasflow -v /u01/DXJ/RNAseq/docker_data:/root/docker_data--privileged=true -it zhxiaokang/rasflow/bin/bash
  3. RASflow项目代码位于docker内的 /root/RASflow

二、进入docker rasflow后conda activate rasflow 并cd /root/RASflow后准备开始运行

图片

（一）按要求整理数据及配置文件(用vim)

  1. 输入、输出数据置于共享卷：主机:容器  /home/tianshu/dxj/docker_data:/home/docker_data

  2. 将原始测序数据上传至主机卷/home/tianshu/dxj/docker_data(容器内即可同步，反之亦然)

    Fastq-103584：包含130个样本(SRR6013474--603)的260个双端原始测序.fastq.gz文件

  3. 配置文件位于容器内/root/RASflow/

    ①数据描述表格：configs/metadata.tsv共三列，sample,group,subject
    ②在configs/config_main.yaml选择参数以控制运行流程，第一次QC yes 第二次no

    ●PROJECT：给项目命名

    ●QC：

    ○第一次先选yes,然后在 output/[PROJECT]/fastqc/report quality control.html查看质控报告

    ○第二次选QC为no：同时根据质控报告决定Trimmed为yes or no；

    ■trimmed后的数据存储在data/output/trim/；trimmed后数据的质控报告在output/[PROJECT]/fastqc after trimming/report quality control after trimming.html

    ●KEY：key file if the data is stored remotely, otherwise leave it tempty

    ●READSPATH：fastq.gz文件存储地

    ●METAFILE: configs/metadata.tsv

    ●END：双端or单端测序

    ●NCORE：首先用getconf _NPROCESSORS_ONLN查看可用核心数，再设定核心数；

    ●REFERENCE：参考为transcriptome or genome，相应地设置参考文件的地址

    ○TRANS: 转录组参考文档所在地

    ○GENOME: 基因组参考文档所在地+ANNOTATION:注释文件所在地+ ATTRIBUTE: gene_id(注释文件中的基因命名法通常为"gene_id",但可能因不是而报错)

    ●ALIGNER: hisat2（当map to genome时方才有用）

    ●COUNTER: featureCounts(默认计数方法)；亦可选htseq-count

    ●DEA: 默认yes

    ●DEATOOL: edgeR【默认】，或 DESeq2

    ●FILTER: 做DEA时，是否将表达值低的基因/转录本过滤掉

    ●GENE_LEVEL：基因水平or转录本做DEA。若GENE_LEVEL为TRUE：

    ○进一步设置ENSEMBL(reference transcriptome)，代表参考转录组是否下载自ENSEMBL

    ■ENSEMBL为TRUE，则通过config/EnsemblDataSet look up table.csv设置EnsemblDataSet

    ■ENSEMBL为FALSE，则在tx2gene中写入两列(tab分割)：transcript ID和对应的gene ID

    ●OUTPUTPATH：存储中间结果(如.bam文件）

    ●FINALOUTPUT：最终结果输出地

    • CONTROL: ["Untreated"]
    • TREAT: ["Dexamethasone"]
    • sample为样本名称(若为双端两个fastq文件，仍然写为SRR6013601即可)，顺序排序
    • group为样本分组信息：表格中的分组信息与config_main.yaml中的信息匹配

    • subject：为配对样本设计所需，通常留空即可，则config_main.yaml中的PAIR信息不会再调用

    （二）python main.py 运行一次判读QC→决定第二次运行是否做数据修剪(trim)

    三、主要的报错是0)metadata.tsv是否正确匹配 1)workflow/*.rules【包括qc,trim,quantify_trans,align】的文件通配符问题：R1/2去掉R和 2)scp改为cp并去掉-i

    ●MissingOutputException in line 35 of /root/RASflow/workflow/trim.rules: ○Missing files after 5 seconds:/root/RASflow/data/output/GSE103584/trim/SRR6013601_*1_val_1.fq.gz/root/RASflow/data/output/GSE103584/trim/SRR6013601_*2_val_2.fq.gz

    ●This might be due to filesystem latency. If that is the case, consider to increase the wait time with --latency-wait.

    写在最后

    单独看笔记可能是不太容易懂的，联系之前的那篇推文分享你的NGS数据分析流程也能发文章哦和RASflow流程原始文献跑一跑可以更好理解。此外，对于不追求一体化、流程化的同学，该流程中各部分使用什么样的软件最合适、输入什么样的文件等问题，也值得学习和探讨。


    收录于合集 #转录组周更
     44个
    上一篇在线/本地获取gmt文件进行GSEA分析下一篇Harvard Chan Bioinformatics Core学习资源介绍
    喜欢此内容的人还喜欢

    微信扫一扫
    关注该公众号