总流程:Fastq->trim->qc->STAR->sort->featureCounts->基因表达矩阵
数据来源:某样本三个生物学重复数据
打开vmware workstation->打开ubuntu操作系统->进入工作目录
激活环境!!!
source ~/miniconda3/bin/activate
cd /mnt/hgfs
cd training_RNA-seq
conda install -c bioconda fastqc=0.12.1
fastqc --version
手动安装软件:
sudo apt-get update
sudo apt-get install build-essential wget
sudo apt install multiqc
使用系统默认版本1.12基本功能相似没有依赖包冲突
trim-galore
github下载压缩包,Linux系统一般下载tar.gz格式
cd /mnt/hgfs/training_RNA-seq/software
tar -xzf TrimGalore-0.6.7.tar.gz
cd TrimGalore-0.6.7
sudo cp trim_galore /usr/local/bin/
sudo chmod +x /usr/local/bin/trim_galore
trim_galore --version
Trim Galore
是一个用于处理高通量测序数据的工具,它实际上依赖于 cutadapt
来进行读长修剪操作。cutadapt
是一个独立的工具,而 Trim Galore
内置了 cutadapt
的一部分功能。
对于 Trim Galore
版本 0.6.7-1 来说,cutadapt
是它的一个依赖项,但它通常会自动处理这个依赖关系。你可以通过以下命令检查 cutadapt
的版本:
cutadapt --version
conda安装依赖
conda install -c bioconda cutadapt
更新conda之后再安装指定版本的cutadapt(只更新不安装cutadapt后续手动安装)
conda update -n base -c conda-forge conda
pip install cutadapt==2.10
手动安装//pip命令未完成我也没安装完pip继续进入software指定cutadapt目录继续操作!!!!
tar -xzvf cutadapt-2.10.tar.gz
cd cutadapt-2.10
pip install .
cutadapt --version
STAR: STAR 2.7.10b
解压进目录
tar -xzf STAR-2.7.10b.tar.gz
cd STAR-2.7.10b/source
编译
make STAR
将编译后的文件移动到系统路径
sudo mv STAR /usr/local/bin/
验证安装
STAR --version
可能的问题
Parameters_openReadsFiles.cpp:25:20: warning: variable ‘imate’ might be clobbered by ‘longjmp’ or ‘vfork’ [-Wclobbered]
编译器提示 Parameters_openReadsFiles.cpp
文件中的变量 imate
可能会被 longjmp
或 vfork
函数的调用覆盖。这是一个编译警告,通常不会导致编译失败,但可能会在特定条件下引发运行时问题。
STAR --help
file STAR
file star应该返回:应该返回类似于 ELF 64-bit LSB executable
的信息,表示文件是一个可执行程序。
可以再进行一次版本安装验证
如果 STAR --version
没有返回任何输出,但没有出现错误,这可能是由于以下几个原因之一:
可能的原因和解决方法:
-
默认情况下不打印版本信息:
STAR --help
2.检查可执行文件是否正常:
file STAR
3.路径问题确保你在运行
STAR --version
时,已经在正确的路径下,或者将STAR
添加到你的环境变量PATH
中,你可以尝试直接调用./STAR --version
(如果你在STAR
可执行文件所在的目录):./STAR --version
4.执行权限问题:确保
STAR
可执行文件具有执行权限。如果没有执行权限,可以使用以下命令赋予权限:chmod +x STAR
hisat2安装
conda install -c bioconda hisat2=2.2.1 hisat2 --version
手动下载featurecounts:
featureCounts
是subread
软件包中的一个工具。虽然featureCounts
是单独的工具,但它确实是subread
包的一部分,所以通常通过安装subread
包来获取featureCounts(注意逐步解决大家的问题)!!!!!
tar -xzvf subread-2.0.6-source.tar.gz cd subread-2.0.6 cd src make -f Makefile.Linux sudo apt-get install build-essential make -f Makefile.Linux clean make -f Makefile.Linux featureCounts -v
bedtools
tar -zxvf bedtools-2.31.0.tar.gz cd bedtools-2.31.0 sudo apt-get update sudo apt-get install libbz2-dev sudo apt-get update sudo apt-get install liblzma-dev sudo make install bedtools --version
Linux和Bash并不完全相同,但它们是密切相关的。
- **Linux** 是一种操作系统内核,它是许多操作系统(如Ubuntu、Debian、CentOS等)的核心部分。Linux提供了操作系统的基础功能,比如硬件控制、文件系统管理、进程管理等。
- **Bash**(Bourne Again Shell)是Linux操作系统中的一种命令行解释器,或者说是一种命令行界面(CLI)。Bash允许用户通过命令行与Linux操作系统进行交互,执行各种任务,例如文件管理、进程控制、脚本编写等。
简单来说,Linux是操作系统的核心,而Bash是与这个系统交互的工具或界面之一。在Linux系统中,Bash是最常用的shell,但也有其他shell,如Zsh、Fish等
关于生物学的部分
之所以想写这部分的原因,是想赋予我的分析生物学意义,毕竟我们所做的一切都是需要有一个研究出口的,于我而言生物信息学重在生物学科,计算机科学在这里只是辅助工具。
上游分析完会得到诸如下图的表达矩阵
这个表达矩阵大有学问:
初步解释:
这个矩阵来自于RNA-seq数据分析,通常是在使用工具如`featureCounts`后生成的基因表达矩阵。以下是各列数据的含义:
1. **Geneid (列B)**:
- 这是基因的标识符(ID),通常使用的是标准化的基因名称或编号。在你的数据中,例如“AT1G01010”是一个基因的标识符。
2. **Chr (列C)**:
- 这是基因所在的染色体编号或名称。例如,“Chr1”表示该基因位于染色体1上。
3. **Start (列D)**:
- 这是基因在染色体上的起始位置,即基因的第一个碱基对在染色体上的位置。
4. **End (列E)**:
- 这是基因在染色体上的终止位置,即基因的最后一个碱基对在染色体上的位置。
5. **Strand (列F)**:
- 这是基因所在的DNA链,通常用“+”表示正义链(sense strand),用“-”表示反义链(antisense strand)。
6. **Length (列G)**:
- 这是基因的长度,通常是基因的起始和终止位置的差值,表示基因占据的碱基对数。
7. **各个样本的表达量 (列H 及其后)**:
- 这些列显示的是不同处理条件下(例如不同的实验重复或者不同的处理组)基因的表达量。每列的标题通常是样本的名称,例如 `cut.Treatment_rep1.sorted.bam` 表示的是一个处理组的第一个重复实验的表达数据。
- 每个单元格中的数字表示的是某个基因在特定样本中的测序读取数,这些数字通常是经过标准化的计数,如TPM(Transcripts Per Million)或FPKM(Fragments Per Kilobase of transcript per Million mapped reads)。
### RNA-seq测到的RNA是什么时期的?
RNA-seq测到的RNA是从样本中提取的mRNA分子。RNA的提取通常是在细胞或组织在特定条件下经过一定时间处理后的某一时点进行的。因此,RNA-seq分析反映的是样本在那个特定时间点或条件下的基因表达情况。至于具体是哪个“时期”的RNA,这取决于实验设计,例如:
- **时间点实验**:如果研究的是发育过程中的基因表达,RNA-seq可能会在多个时间点进行,以观察基因表达随时间的变化。
- **处理条件**:如果研究的是药物或环境条件的影响,RNA可能是在施加处理一段时间后提取的,以反映处理后的基因表达。
总之,RNA-seq测到的RNA对应的是你在实验设计中选择的那个时间点或条件下细胞内活跃的mRNA,代表了那个特定时期的基因表达谱。
更详细理解
基因的起始位置和终止位置
在基因组中,基因的起始位置和终止位置通常并不只是一个单一的数字,因为基因可能由多个**外显子**(exons)组成,每个外显子在染色体上都有各自的起始和终止位置。这些外显子通过**内含子**(introns)分隔开。因此,对于一个基因来说,起始位置和终止位置的定义可以分为以下几种情况:
1. **基因水平的起始和终止位置**:
- **起始位置**:这是整个基因的第一个外显子的第一个碱基在染色体上的位置。
- **终止位置**:这是整个基因的最后一个外显子的最后一个碱基在染色体上的位置。
举个例子,如果一个基因由三个外显子组成,且它们在染色体上的位置分别是:
- 外显子1:起始位置 100, 终止位置 200
- 外显子2:起始位置 500, 终止位置 600
- 外显子3:起始位置 800, 终止位置 900
那么这个基因的起始位置就是 100,终止位置就是 900。
2. **外显子水平的起始和终止位置**:
- 每个外显子都有自己的起始和终止位置。如果你在表格中看到多个位置,例如`3631:3913, 3991:4276`,这些通常表示这个基因有多个外显子,每个外显子的起始和终止位置都被列出。
### 基因长度的计算
基因长度通常有两种计算方法:
1. **基因组水平的基因长度**:
- 这是从第一个外显子的起始位置到最后一个外显子的终止位置的距离。例如,若基因起始于 100,结束于 900,则长度为 900 - 100 + 1 = 801 个碱基对。
- 注意,这种计算包括了内含子,也就是计算的是整个基因在基因组中占据的区间。
2. **外显子总和长度**:
- 这是将所有外显子的长度相加的总和,不包括内含子。例如,如果一个基因有三个外显子,其长度分别为 100、101 和 101 碱基对,那么基因的总长度为 100 + 101 + 101 = 302 个碱基对。
- 这种计算方法通常用来表示基因编码区域的总长度,也就是基因实际能翻译成蛋白质的区域长度。
数据中的示例
在上图的矩阵中,长度一栏代表的就是上面提到的基因长度。比如,对于基因 AT1G01010(在 Chr1 上),其外显子起始和终止位置是`3631:3913, 3991:4276`,这些位置的长度就是每个外显子长度的总和。
假设第一个外显子的长度是 `3913 - 3631 + 1 = 283` 碱基对,第二个外显子的长度是 `4276 - 3991 + 1 = 286` 碱基对,那么整个基因的外显子总和长度就是 `283 + 286 = 569` 碱基对。如果这个基因的长度被记录为更大值,可能表示的是基因的整个区间长度,而不是外显子总和长度。
### 总结
- **起始和终止位置**可以是多个,表示基因中多个外显子的范围。
- **基因长度**可以表示整个基因组中的区间长度也可以理解为外显子总和长度。