bulk-RNA seq测序数据分析流程

假如有bulk-RNA测序的数据:TH1,TH2,TH3三个重复(实验组),TW1,TW2,TW3三个重复(对照组)

  • 准备工作
需要安装的软件(如FastQC、Trimmomatic、HISAT2、StringTie、samtools)
conda install -c bioconda fastqc
conda install trimmomatic
conda install -c bioconda stringtie 
conda install -c bioconda hisat2
conda install samtools
conda install -c bioconda picard

1. 质量控制

首先,需要对原始的测序数据(FASTQ格式)进行质量控制。通常使用FastQC进行初步的质量检查,然后使用Trimmomatic或者其他工具进行质量修剪。

# 质量检查
fastqc *.fq.gz# 质量修剪
for sample in TH1 TH2 TH3 TW1 TW2 TW3
dotrimmomatic PE ${sample}.R1.fq.gz ${sample}.R2.fq.gz \${sample}.R1.trim.fq.gz ${sample}.R1.untrim.fq.gz \${sample}.R2.trim.fq.gz ${sample}.R2.untrim.fq.gz \SLIDINGWINDOW:4:20 MINLEN:25
done

2.读段比对

使用比对工具(如HISAT2, STAR等)将质量修剪后的读段比对到参考基因组。

  • 使用HISAT2进行读段比对
# 假设参考基因组索引为genome_index
for sample in TH1 TH2 TH3 TW1 TW2 TW3
dohisat2 -x genome_index -1 ${sample}.R1.trim.fq.gz -2 ${sample}.R2.trim.fq.gz \-S ${sample}.sam
done
  • 使用STAR进行读段比对
  1. 下载并构建STAR索引

在使用STAR进行比对之前,需要首先构建或下载参考基因组的索引。这是一个一次性的过程。假设你已经有了参考基因组的fasta文件和基因注释文件(GTF格式),使用以下命令构建索引:

STAR --runThreadN NumberOfThreads \
--runMode genomeGenerate \
--genomeDir /path/to/genomeDir \
--genomeFastaFiles /path/to/genome/fasta \
--sjdbGTFfile /path/to/annotations.gtf \
--sjdbOverhang ReadLength-1

其中,NumberOfThreads是使用的线程数,/path/to/genomeDir是存储索引的目录,/path/to/genome/fasta和/path/to/annotations.gtf分别是参考基因组的fasta文件和基因注释文件的路径。ReadLength-1是读段长度减去1,这是一个重要的参数,需要根据数据来设定。
2. 读段比对

for sample in TH1 TH2 TH3 TW1 TW2 TW3
doSTAR --genomeDir /path/to/genomeDir \--readFilesIn ${sample}.R1.trim.fq.gz ${sample}.R2.trim.fq.gz \--runThreadN 10 \--outFileNamePrefix ${sample} \--readFilesCommand zcat \--outSAMtype BAM SortedByCoordinate
done

在这里,/path/to/genomeDir是你的索引目录,NumberOfThreads是使用的线程数, s a m p l e . R 1. t r i m . f q . g z 和 {sample}.R1.trim.fq.gz和 sample.R1.trim.fq.gz{sample}.R2.trim.fq.gz分别是质量修剪后的读段文件。–outSAMtype BAM SortedByCoordinate会输出排序后的BAM文件,这对后续分析非常有用。

  • HISAT2和STAR的区别如下:

    1. 设计和算法:
      HISAT2:使用了一种称为分层索引的方法。它首先将参考基因组划分为较小的片段,并为这些片段构建全局和局部的索引。这种方法使得HISAT2在处理基因组比对时特别有效,尤其是在考虑剪接位点和变异位点时。
      STAR:使用了一种称为基于稀疏后缀数组的方法。STAR的这种方法允许它在一个单独的步骤中完成比对和剪接位点的检测,这使得它在处理大型数据集时非常快速和内存高效。

    2. 性能和资源消耗:
      HISAT2:通常比STAR使用更少的内存,但在处理时间上可能稍长。它特别适合于资源受限的计算环境。
      STAR:虽然需要较多的内存,但在执行速度上通常比HISAT2更快,特别是在处理包含大量剪接事件的复杂转录组时。

    3. 准确性和灵敏度:
      HISAT2:在处理较长的读段和考虑基因组变异(如SNPs和小的indels)时表现出较高的准确性。
      STAR:在检测剪接事件和转录本多样性方面表现优异,尤其是在处理短读段数据时。

    4. 适用性: HISAT2:由于其对于基因组变异的考虑,HISAT2在进行种群遗传学研究或者需要高精度比对的场景中表现更好。 STAR:在进行转录组学研究,尤其是需要快速处理大规模数据集和复杂剪接事件的情况下更为适用。

    5. 易用性和功能: 两者都具备易用性和丰富的功能集,包括对多线程的支持、与常用转录组分析流程的兼容性等。

综上,选择HISAT2还是STAR取决于你的具体需求,比如数据的类型、可用资源、分析的焦点等。在某些情况下,研究者可能会根据不同阶段的需求选择不同的工具。例如,在初步筛选或快速分析时使用STAR,而在需要更细致分析基因组变异或长读段数据时使用HISAT2。

3.标记和去除PCR重复

使用Picard或samtools来标记和去除重复的读段。
对于使用samtools,流程如下:

for sample in TH1 TH2 TH3 TW1 TW2 TW3
dosamtools view -bS ${sample}.sam | \samtools sort -o ${sample}.sorted.bamsamtools markdup ${sample}.sorted.bam ${sample}.dedup.bam
done

对于使用Picard, 流程如下:

for sample in TH1 TH2 TH3 TW1 TW2 TW3
dosamtools view -bS ${sample}.sam | \samtools sort -o ${sample}.sorted.bampicard MarkDuplicates I=${sample}.sorted.bam \O=${sample}.dedup.bam M=${sample}.marked_dup_metrics.txt
done

4.转录组组装和表达量估计

使用工具如StringTie或Cufflinks对SAM/BAM文件进行转录组组装,生成转录本。


for sample in TH1 TH2 TH3 TW1 TW2 TW3
dostringtie -e -B -G reference.gtf -o ${sample}.gtf ${sample}.dedup.bam
done

-e 选项告诉StringTie只估计已知转录本(在提供的GTF文件中)的表达量。
-B 选项会生成每个基因和转录本的表达量计数文件,这些文件可用于后续差异表达分析。
-G reference.gtf 是参考基因组的注释文件,这里需要你提供正确的路径和文件名。

5.差异表达分析

最后,使用DESeq2或edgeR等工具进行差异表达分析,以识别实验组和对照组之间差异表达的基因。

library(DESeq2)# 假设countData为一个矩阵,行为基因,列为样本,包含基因的计数
# colData为DataFrame,包含样本信息
# 这里需要根据实际情况准备countData和colDatadds <- DESeqDataSetFromMatrix(countData = countData,colData = colData,design = ~ condition)dds <- DESeq(dds)
res <- results(dds)# 保存结果
write.csv(as.data.frame(res), file="DESeq2_results.csv")

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/601280.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

adb shell getevent/sendevent

#### adb shell getevent 获取点击事件 100135925:/ # getevent add device 1: /dev/input/event2name: "mtk-tpd" /dev/input/event2: 0001 014a 00000001 /dev/input/event2: 0003 0039 00000088 /dev/input/event2: 0003 0035 00000072 /dev/input/event2: 00…

【Pytorch】学习记录分享10——TextCNN用于文本分类处理

【Pytorch】学习记录分享10——PyTorchTextCNN用于文本分类处理 1. TextCNN用于文本分类2. 代码实现 1. TextCNN用于文本分类 具体流程&#xff1a; 2. 代码实现 # coding: UTF-8 import torch import torch.nn as nn import torch.nn.functional as F import numpy as np…

14、接口

接口 ​ 接口interface&#xff0c;是一组行为规范的集合&#xff0c;就是定义一组未实现的函数声明。谁使用接口就是参照接口的方法定义实现它们。 type 接口名 interface {方法1 (参数列表1) 返回值列表1方法2 (参数列表2) 返回值列表2... }接口命名习惯在接口名后面加上er…

【机器学习:欧氏距离 】机器学习中欧氏距离的理解和应用

【机器学习&#xff1a;欧氏距离 】机器学习中欧氏距离的理解和应用 距离公式二维更高的维度点以外的物体属性欧几里得距离的平方概括历史 在数学中&#xff0c;欧氏距离’是指欧氏空间中任意两点之间的直线距离。这种距离可以通过应用勾股定理来计算&#xff0c;利用两点的笛卡…

如何停止一个运行中的Docker容器

要停止一个运行中的Docker容器&#xff0c;你可以使用以下命令&#xff1a; docker stop <容器ID或容器名> 将 <容器ID或容器名> 替换为你要停止的具体容器的标识符或名称。你可以使用以下命令查看正在运行的容器&#xff1a;docker ps 这将列出所有正在运行的…

Linux内核(2)-Makefile详解,必须要掌握的编译参数

1.版本号 VERSION 4 PATCHLEVEL 1 SUBLEVEL 152.MAKEFLAGS变量 MAKEFLAGS -rR --include-dir$(CURDIR) 包含当前目录及所有递归子目录 3.make V1编译输出 make V1 输出编译完整命令 ifeq ("$(origin V)", "command line")KBUILD_VERBOSE $(V) en…

再检查下这些测试思维面试题你都会了么?

创建坐席组的功能模块&#xff0c;如何进行测试用例设计&#xff1f; 解答&#xff1a; 功能测试&#xff0c;使用等价类划分法去分析创建坐席的每个输入项的有效及无效类&#xff0c;同步考虑边界值去设计对应的测试用例&#xff1a; 先进行冒烟测试&#xff0c;正常创建坐席…

操作系统期末复习知识点

目录 一.概论 1.操作系统的介绍 2.特性 3.主要功能 4.作用 二.进程的描述与控制 1.进程的定义 2.特性 3.进程的创建步骤 4.基本状态转化 5.PCB的作用 6.进程与线程的比较 三.进程同步 1.同步的概念&#xff08;挺重要的&#xff09; 2.临界区 3.管程和进程的区…

62.网游逆向分析与插件开发-游戏增加自动化助手接口-游戏公告类的C++还原

内容来源于&#xff1a;易道云信息技术研究院VIP课 上一个内容&#xff1a;游戏红字公告功能的逆向分析-CSDN博客 码云地址&#xff08;master分支&#xff09;&#xff1a;https://gitee.com/dye_your_fingers/sro_-ex.git 码云版本号&#xff1a;0888e34878d9e7dd0acd08ef…

Redis第四讲——Redis的数据库结构、删除策略及淘汰策略

一、redis中的数据库 redis服务器将所有数据库都保存在服务器状态redis.h/redisServer结构的db数组中。db数组的每项都是一个redis.h/redisDb结构&#xff0c;而每个redisDb结构就代表一个数据库。在初始化服务器时&#xff0c;程序会根据服务器状态的dbnum属性来决定应该创建多…

C++ 操作重载与类型转换

文章目录 基本概念为什么使用操作重载&#xff1f;注意事项&#xff1a; 输入和输出运算符重载输出运算符 <<重载输入运算符 >> 算术和关系运算符算术运算符相等运算符关系运算符 赋值运算符基本原则&#xff1a;示例&#xff1a;注意事项&#xff1a; 下标运算符实…

outlook邮件群发单显技巧?群发怎么单显?

outlook邮件群发单显如何设置&#xff1f;QQ邮箱怎么群发单显&#xff1f; 在群发邮件时&#xff0c;如何让每个收件人只看到自己的名字&#xff0c;而不是其他人的名字&#xff0c;这就涉及到所谓的“单显”技巧。下面蜂邮EDM就为大家揭秘Outlook邮件群发单显的奥秘。 outlo…

实战干货:用 Python 批量下载百度图片!

为了做一个图像分类的小项目&#xff0c;需要制作自己的数据集。要想制作数据集&#xff0c;就得从网上下载大量的图片&#xff0c;再统一处理。 这时&#xff0c;一张张的保存下载&#xff0c;就显得很繁琐。那么&#xff0c;有没有一种方法可以把搜索到的图片直接下载到本地电…

Spark 运行架构

Spark 框架的核心是一个计算引擎&#xff0c;整体来说&#xff0c;它采用了标准 master-slave 的结构。 如下图所示&#xff0c;它展示了一个 Spark 执行时的基本结构。图形中的 Driver 表示 master&#xff0c; 负责管理整个集群中的作业任务调度。图形中的 Executor 则是 sla…

Swift爬虫使用代理IP采集唯品会商品详情

目录 一、准备工作 二、代理IP的选择与使用 三、使用Swift编写唯品会商品爬虫 四、数据解析与处理 五、注意事项与优化建议 六、总结 一、准备工作 在开始编写爬虫之前&#xff0c;需要准备一些工具和库&#xff0c;以确保数据抓取的顺利进行。以下是所需的工具和库&…

速盾高防ip:专业防御ddos

速盾高防IP是速盾网络为企业提供的专业DDoS攻击防御解决方案之一。作为一种先进的网络安全服务&#xff0c;速盾高防IP致力于保护客户的网络资源免受分布式拒绝服务&#xff08;DDoS&#xff09;攻击的威胁。以下是速盾高防IP的一些关键特点和优势&#xff1a; 实时攻击监测&am…

Mac 软件出现「意外退出」及「打不开」解决方法

Mac 软件出现「意外退出」及「打不开」解决方法 软件出现意外退出及软件损坏的情况&#xff0c;这是因为苹果删除了TNT的证书&#xff0c;所以大部分TNT破解的Mac软件会出现无法打开&#xff0c;提示意外退出。 终端需先安装Xcode或Apple命令行工具 如未装Xcode可以使用下列命…

【算法设计与分析】分治-时间复杂度计算

目录 主定理 Master Theorem分治算法运行时间的递归表示主定理的简化形式 主定理的一般形式 递归树 Recursion Tree递归树的简单结论 主定理 Master Theorem 分治算法运行时间的递归表示 将原问题分解成 a 个子问题递归求解&#xff0c;每个子问题的规模是原问题的 1/b。同时子…

go-cqhttp作者停止维护——替代品OpenShamrock的使用方法

目录 前言 解决办法 配置要求 实操 刷入面具 安装lsp框架 安装OpenShamrock和QQ 注意 大功告成 前言 由于QQ官方针对协议库的围追堵截&#xff0c;go-cqhttp已经无力维护下去了 原文连接 QQ Bot的未来以及迁移建议 Issue #2471 Mrs4s/go-cqhttp (github.com)https…

libcurl的get、post的使用

demo使用的是curl-8.3.0.tar.gz&#xff0c;其它版本也可以&#xff0c;安装教程可以去网上搜 #include <stdio.h> #include <stdlib.h> #include <string.h> #include <errno.h> #include <assert.h> /* somewhat unix-specific */ #include &…