bs自用教程

sepasure

已于 2024-07-29 22:35:40 修改

阅读量414

点赞数 11

文章标签： python r语言

于 2024-07-24 17:53:38 首次发布

本文链接：https://blog.csdn.net/qq_45144772/article/details/140574370

版权

前言

参考https://github.com/erilu/bulk-rnaseq-analysis?tab=readme-ov-file#concatenating-fastq-files
以下示例：我们的目标是找到差异表达的基因，以响应 LNCaP 和 PC3 细胞系的缺氧。因此，我们将选择两种细胞系的对照样品（Empty_Vector LNCaP 和 siCtrl 的对照样品对于PC3），在常氧或缺氧条件下。具体下表概述了我们需要下载的示例：
NCBI以序列读取存档的形式存储测序数据（SRA）文件。每个样本都与一组 SRA 相关联入藏号，如上所示。首先，我们需要下载SRA 针对每个示例运行。然后，我们将使用 SRA 文件生成 FASTQ 文件。FASTQ文件是批量所需的数据格式 RNA测序分析。

要将 SRA 文件下载到我们的机器上，我们将使用 NCBI 的 SRA 工具箱。

在这里，我使用Linux下载，由于超算账户限制，我在Windows下载好，上传到超算用户文件夹下。具体下载请查看https://github.com/ncbi/sra-tools/wiki/01.-Downloading-SRA-Toolkit

查找数据库

有关在 GEO 中导航并下载 SRA Runs，可以查看https://erilu.github.io/python-fastq-downloader/
已从生物公司获取数据，直接跳转使用 STAR 映射读取

FASTQ的获取

windows上直接fastq下载https://www.ncbi.nlm.nih.gov/Traces/index.html?view=run_browser&display=metadata
linux下载：（超算）

单个下载sra：

conda activate bulkseq
cd /public/home/kout/sra# 提前cd到想要下载的文件夹路径
prefetch SRR7179504

运行结果

(bulkseq) [kout@tc6000 sra]$ 2020-02-14T22:15:36 prefetch.2.8.2: 1) Downloading 'SRR7179504'...
(bulkseq) [kout@tc6000 sra]$ 2020-02-14T22:15:36 prefetch.2.8.2:  Downloading via https...
(bulkseq) [kout@tc6000 sra]$ 2020-02-14T22:18:39 prefetch.2.8.2: 1) 'SRR7179504' was downloaded successfully
(bulkseq) [kout@tc6000 sra]$ 2020-02-14T22:18:39 prefetch.2.8.2: 'SRR7179504' has 0 unresolved dependencies

创建 SRR7179504 的 FASTQ 文件

fastq-dump --outdir fastq --gzip --skip-technical  --readids --read-filter pass --dumpbase --split-3 --clip /public/home/kout/sra/SRR7179504

结果如下

Rejected 13548432 READS because of filtering out non-biological READS
Read 13548432 spots for /public/home/kout/sra/SRR7179504
Written 13548432 spots for /public/home/kout/sra/SRR7179504

批量下载sra：

python fastq_download.py

以下是fastq_download.py

import subprocess

# 提前cd到想要下载的文件夹路径
sra_numbers = [
    "SRR7179504", "SRR7179505", "SRR7179506", "SRR7179507",
    "SRR7179508", "SRR7179509", "SRR7179510", "SRR7179511",
    "SRR7179520", "SRR7179521", "SRR7179522", "SRR7179523",
    "SRR7179524", "SRR7179525", "SRR7179526", "SRR7179527",
    "SRR7179536", "SRR7179537", "SRR7179540","SRR7179541"
    ]

for sra_id in sra_numbers:
    print ("Currently downloading: " + sra_id)
    prefetch = "prefetch " + sra_id
    print ("The command used was: " + prefetch)
    subprocess.call(prefetch, shell=True)

# this will extract the .sra files from above into a folder named 'fastq'
for sra_id in sra_numbers:
    print ("Generating fastq for: " + sra_id)
    fastq_dump = "fastq-dump --outdir fastq --gzip --skip-technical  --readids --read-filter pass --dumpbase --split-3 --clip /public/home/kout/sra/" + sra_id + "/" + sra_id + ".sra"
    print ("The command used was: " + fastq_dump)
    subprocess.call(fastq_dump, shell=True)

连接 FASTQ 文件
每个 LNCaP 样本都与四次 SRA 运行相关联，对于每个样本，我们应该将使用命令将四个文件合并到一个 FASTQ 文件中。下面，我对每个 LNCaP 样本执行串联：

cat SRR7179504_pass.fastq.gz SRR7179505_pass.fastq.gz SRR7179506_pass.fastq.gz SRR7179507_pass.fastq.gz  > LNCAP_Normoxia_S1.fastq.gz
cat SRR7179508_pass.fastq.gz SRR7179509_pass.fastq.gz SRR7179510_pass.fastq.gz SRR7179511_pass.fastq.gz  > LNCAP_Normoxia_S2.fastq.gz
cat SRR7179520_pass.fastq.gz SRR7179521_pass.fastq.gz SRR7179522_pass.fastq.gz SRR7179523_pass.fastq.gz  > LNCAP_Hypoxia_S1.fastq.gz
cat SRR7179524_pass.fastq.gz SRR7179525_pass.fastq.gz SRR7179526_pass.fastq.gz SRR7179527_pass.fastq.gz  > LNCAP_Hypoxia_S2.fastq.gz

相比之下，每个 PC3 样本只有一个 FASTQ 文件。我们可以将它们从其 SRR 标识符重命名为其真实样本使用以下名称：

mv SRR7179536_pass.fastq.gz PC3_Normoxia_S1.fastq.gz
mv SRR7179537_pass.fastq.gz PC3_Normoxia_S2.fastq.gz
mv SRR7179540_pass.fastq.gz PC3_Hypoxia_S1.fastq.gz
mv SRR7179541_pass.fastq.gz PC3_Hypoxia_S2.fastq.gz

单个 SRA 可以将它们全部删除。
现在，该文件夹应包含总共 8 个 FASTQ 文件：4 个用于 LNCaP，4 个用于 PC3。我们已准备好开始调整这些内容 FASTQ 文件复制到参考基因组！

使用 STAR 映射读取

映射读取背后的概念
FASTQ 文件中的每个读段都必须“映射”到参考基因组。 “定位”可以描述为在参考基因组中找到位置。

为了进行批量RNA测序，“文库”必须首先从目标细胞的RNA制备。这个库是排序以创建 FASTQ 文件。要创建文库，mRNA 来自首先捕获细胞并逆转录成cDNA。这就是碎片化成小块，并将测序适配器连接到每个片段的末端以创建库。音序器将“读取” 这些片段并将序列存储在 FASTQ 文件中。

我们可以打印出其中一个“reads”来显示序列的样子，存储在 FASTQ 文件中的 4 行块中。

zcat LNCAP_Hypoxia_S1.fastq.gz | head -4

@SRR7179520.1.1 1 length=76
GTGAANATAGGCCTTAGAGCACTTGANGTGNTAGNGCANGTNGNNCCGGAACGNNNNNNNNAGGTNGNNNGNGTTG
+SRR7179520.1.1 1 length=76
AAAAA#EEEEEEEEEEEEEEEEEEEE#EEE#EEE#EEE#EE#E##EEEEEEEE########EEEE#E###E#EAEA

读取的每一行都包含一组不同的信息：

第 1 行包含带有可选说明的序列标识符。此行必须以“@”开头。
第 2 行包含读取的实际顺序。
第 3行必须以“+”开头，并且通常包含相同的序列标识符为第 1 行。不需要具有标识符，并且有时，此行仅由一个“+”组成。
第 4 行包含ine 2 中每个基本名称的质量值。例如，第一个碱基“G”的质量值为“A”。这是基本上是对机器的信心评估分配了正确的基地。

使用每个片段的序列（第 3 行），我们可以使用软件来弄清楚该片段最初来自基因组的哪一部分，基于序列相似性。这告诉我们哪个基因正在被转录以创建原始 mRNA，该 mRNA 被捕获以制造图书馆。计算读取映射到特定基因的次数为我们提供了有关基因“高”或“低”的信息表示。

建立基因组索引
映射读取的一个很好的工具是 STAR。https://github.com/alexdobin/STAR

以下为我的记录;

wget https://github.com/alexdobin/STAR/archive/2.7.0e.tar.gz
tar -xzf 2.7.0e.tar.gz
cd STAR-2.7.0e
cd /public/home/kout/sra/fastq/STAR-2.7.0e/source
make STAR

接下来要将STAR写入环境变量中
如果不写环境变量会报错找不到STAR，参考https://github.com/alexdobin/STAR/issues/1001

vim ~/.bashrc

在最后一行插入

export PATH=$PATH:/public/home/kout/sra/fastq/STAR-2.7.0e/source

使文件配置立即生效

source ~/.bashrc

注意生效后会自动切换base环境

conda activate bulkseq

我将使用版本 2.7.0e. 安装 STAR 后，我们需要首先建立一个基因组索引。由于我们正在分析来自人类细胞的 RNA-seq 数据，因此我们将构建人类参考基因组。

为了建立基因组索引，我们必须在模式下运行 STAR，并为其提供人类参考基因组序列（FASTA）和注释（GTF）文件。

wget ftp://ftp.ensembl.org/pub/release-99/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
wget ftp://ftp.ensembl.org/pub/release-99/gtf/homo_sapiens/Homo_sapiens.GRCh38.99.gtf.gz
gunzip *.gz

接下来制作基因组索引文件，指定保存 FASTA 和 GTF 文件的目录。
运行这句需要大量的内存，我在超算中好几次都是在sorting Suffix Array chunks and saving them to disk...后进程被killed，尝试slurm提交并申请64G内存可以跑通。以下为我的slurm脚本：

#!/bin/bash
#SBATCH -J bulkseq  
#SBATCH -o bulkseq.out   
#SBATCH -e bulkseq.err    
#SBATCH -N 1                         
#SBATCH -n 1                 
#SBATCH --time=4:00:00                    
#SBATCH -p ojfat  
#SBATCH --mem=64G
         
STAR --runThreadN 16 \
     --runMode genomeGenerate \
     --genomeDir //public/home/kout/sra/fastq/STAR-2.7.0e \
     --genomeFastaFiles //public/home/kout/sra/fastq/STAR-2.7.0e/Homo_sapiens.GRCh38.dna.primary_assembly.fa \
     --sjdbGTFfile //public/home/kout/sra/fastq/STAR-2.7.0e/Homo_sapiens.GRCh38.99.gtf \
     --genomeSAsparseD 2

输出如下所示：

Jul 27 21:55:23 ..... started STAR run
Jul 27 21:55:23 ... starting to generate Genome files
Jul 27 21:56:19 ... starting to sort Suffix Array. This may take a long time...
Jul 27 21:56:30 ... sorting Suffix Array chunks and saving them to disk...
Jul 27 22:50:54 ... loading chunks from disk, packing SA...
Jul 27 22:51:35 ... finished generating suffix array
Jul 27 22:51:35 ... generating Suffix Array index
Jul 27 22:53:39 ... completed Suffix Array index
Jul 27 22:53:39 ..... processing annotations GTF
Jul 27 22:53:53 ..... inserting junctions into the genome indices
Jul 27 22:59:47 ... writing Genome to disk ...
Jul 27 22:59:49 ... writing Suffix Array to disk ...
Jul 27 22:59:54 ... writing SAindex to disk
Jul 27 22:59:55 ..... finished successfully

运行命令后，您应该会看到一堆基因组索引文件在指定的目录中。我们现在可以开始映射了，从 FASTQ 文件读取到这个新创建的参考基因组。

映射读取

由于我们有 8 个样本要运行，并且每个样本 sample 需要很长时间才能运行，我们可以使用 align_STAR.py 脚本来运行目录中每个 FASTQ 的命令：

python align_STAR.py

以下为python align_STAR.py：

import subprocess
import os

output_directory = "/public/home/linql/STAR-2.7.0e/LNCAP"
genome_directory = "/data/genomes/h38/STAR/"
fastq_directory = "/data/analysis/hypoxia/fastq/"

os.mkdir(output_directory)
for fastq in os.listdir(fastq_directory):
    # only process files that end in fastq.gz
    if fastq.endswith('.fastq.gz'):
        prefix=fastq.strip(".fastq.gz") + "_output"
        # make an output folder for the current fastq file
        os.mkdir(output_directory + prefix)
        print ("Currently mapping: " + fastq)
        # run STAR on the current fastq file
        subprocess.call("STAR --runThreadN 16 --genomeDir " + genome_directory + " --readFilesCommand zcat --outFilterType BySJout --outFilterMismatchNoverLmax 0.04 --outFilterMismatchNmax 999 --alignSJDBoverhangMin 1 --alignSJoverhangMin 8 --outFilterMultimapNmax 20 --alignIntronMin 20 --alignIntronMax 1000000 --alignMatesGapMax 1000000 --readFilesIn "+ fastq_directory + fastq + " --clip3pAdapterSeq GATCGGAAGAGCACACGTCTGAACTCCAGTCAC --outSAMtype BAM SortedByCoordinate --quantMode GeneCounts --outFileNamePrefix " + output_directory + prefix + "/", shell=True)

以下是我将python改写成slurm脚本：

#!/bin/bash
#SBATCH --job-name=STAR_mapping
#SBATCH --output=STAR_mapping.out
#SBATCH --error=STAR_mapping.err
#SBATCH --nodes=1                         
#SBATCH --ntasks=1  
#SBATCH --cpus-per-task=16
#SBATCH --time=2:00:00
#SBATCH --mem=64G
#SBATCH --partition=ojfat  

output_directory="/public/home/linql/STAR-2.7.0e/LNCAP/"
genome_directory="/public/home/linql/STAR-2.7.0e"  # 基因组索引文件
fastq_directory="/public/home/linql/STAR-2.7.0e/LNCAP/FASTAQ"

mkdir -p $output_directory

for fastq in $(ls $fastq_directory); do
    # 仅处理以 fastq.gz 结尾的文件
    if [[ $fastq == *.fastq.gz ]]; then
        prefix=$(basename $fastq .fastq.gz)_output
        mkdir -p ${output_directory}${prefix}
        echo "Currently mapping: $fastq"
        STAR --runThreadN 16 \
             --genomeDir $genome_directory \
             --readFilesCommand zcat \
             --outFilterType BySJout \
             --outFilterMismatchNoverLmax 0.04 \
             --outFilterMismatchNmax 999 \
             --alignSJDBoverhangMin 1 \
             --alignSJoverhangMin 8 \
             --outFilterMultimapNmax 20 \
             --alignIntronMin 20 \
             --alignIntronMax 1000000 \
             --alignMatesGapMax 1000000 \
             --readFilesIn ${fastq_directory}/${fastq} \
             --clip3pAdapterSeq GATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
             --outSAMtype BAM SortedByCoordinate \
             --quantMode GeneCounts \
             --outFileNamePrefix ${output_directory}${prefix}/
    fi
done

运行结果

Currently mapping: PC3_Normoxia_S2.fastq.gz
Jul 28 17:34:44 ..... started STAR run
Jul 28 17:34:45 ..... loading genome
Jul 28 17:34:52 ..... started mapping
Jul 28 17:37:53 ..... finished mapping
Jul 28 17:37:54 ..... started sorting BAM
Jul 28 17:39:30 ..... finished successfully

注意的参数：

--genomeDir，表示包含 Genome Index文件
--readFilesIn，指定示例的 FASTQ 文件位置 LNCAP_Normoxia_S1_R1_001.fastq.gz
--outFileNamePrefix，它指定了我们要给 STAR输出

其余参数为STAR推荐的默认设置，这通常就足够了。

了解 ReadsPerGene.out.tab
在目录中，您将找到一个单独的文件夹您映射的每个 FASTQ 文件的结果。在每个文件夹中，您应看到以下文件：
在这里插入图片描述

这文件包含映射到中每个基因的读取次数转录组。使用以下命令快速浏览一下文件即可显示文件的结构。在开始时有一些QC指标，转录组中每个基因后跟的文件（N_unmapped等）。

cat ReadsPerGene.out.tab | head -n 10

结果：

N_unmapped  2461926 2461926 2461926
N_multimapping  4255262 4255262 4255262
N_noFeature 4048236 41792820    4475290
N_ambiguous 3081395 50727   1691846
ENSG00000223972 0   0   0
ENSG00000227232 3   0   3
ENSG00000278267 13  0   13
ENSG00000243485 0   0   0
ENSG00000284332 0   0   0
ENSG00000237613 0   0   0

N_unmapped: 这是未能成功映射到参考基因组或转录组的读数数量。
N_multimapping: 这是多重映射的读数数量。多重映射意味着一个读数可以映射到多个位置。
N_noFeature: 这是没有映射到任何特征（如基因、外显子等）的读数数量。
N_ambiguous: 这是由于映射位置不明确而无法分配到特定特征的读数数量。
每行表示一个基因（用基因的ENSG编号表示）。
每个基因后面有三个数值，可能表示在三个不同样本或条件下的读数计数。例如：
ENSG00000223972 在所有三个样本中的表达量都是 0。
ENSG00000227232 在第一个和第三个样本中的表达量分别是 3 和 3，在第二个样本中的表达量是 0。

第一列是基因 ID。
第 2、3 和 4 列对应于映射到每个基因的读取数：

Unstranded：未区分正负链的读取数。这列的数值是从两个链上读取数的总和。
Stranded：链特异性正链的读取数。
ReverseStranded：链特异性负链的读取数。

每列对应于读取映射的“搁浅性”。组成一个RNA-seq文库 DNA，它有两条链：“感觉”链和“反义”链链。RNA-seq文库可以制备为“非链”或 “搁浅”。您为下游分析选择的色谱柱通常是由用于制备的文库套件的搁浅性决定样品。

如果使用“非链式”文库制备试剂盒，则第 2 列应为选择。对于未搁浅的库，我们不知道是哪一链原始 mRNA 由两者产生，并且读数映射到两者意义和反意义链。
如果使用了“搁浅”文库制备试剂盒，我们必须从中选择第 3 或第 4 列。mRNA起源于哪条链的信息 from 将保留。第 3 列对应于第一个读取链对齐，第 4 列对应于对齐的第二个读取链。

为了确定我们样品的搁浅度，我们需要找出答案所使用的确切文库制备试剂盒。样本信息 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM3145522 指出，TruSeq® 链式 mRNA 文库制备试剂盒（RS-122-2101， Illumina）用于制备这些样品。由于这是一个搁浅的文库制备试剂盒，我们必须从第 3 列或第 4 列中进行选择。查看上面的输出，我们观察到第 4 列包含读取，而列 3 没有。我们将继续使用每个样本的第 4 列用于下游分析。

按样本选择读取计数

对于我们的每个样本，我们需要从中选择第 4 列，并将它们放入一个数据文件中以供下游使用分析。我们可以编写一个 Python 脚本，而不是手动执行此操作要循环遍历其中的每个文件夹，请获取并将数据写入单个输出文件。该脚本通过存储字典中每个样本的每个基因计数，以及将聚合数据写入名为的文件中。

ReadsPerGene.out.tabSTAR_output/ReadsPerGene.out.tabcombine_STAR_output.pyraw_counts.csv

以下为combine_STAR_output.py


import os
import csv

# modify desired_column based on strandedness of library kit
desired_column = 3
output_directory = "/public/home/linql/STAR-2.7.0e/LNCAP/"
sample_dict = {}
sample_names = []

# loop through each folder and store data from ReadsPerGene.out.tab in sample_dict
for folder in os.listdir(output_directory):
    if folder.endswith("_output"):
        file_name=folder.strip('_output')
        sample_names.append(file_name)
        gene_dict = {}
        with open(output_directory + folder + "/ReadsPerGene.out.tab") as tabfile:
            print("Column " + str(desired_column) + " of " + file_name + " stored")
            reader = csv.reader(tabfile, delimiter = "\t")
            for row in reader:
                gene_dict[row[0]]=row[desired_column]
        # store gene_dict for the current file in sample_dict
        sample_dict[file_name] = gene_dict

# write sample_dict to output files
# the qc metrics are stored in qc.csv, and the gene counts are stored in raw_counts.csv
with open("raw_counts.csv", "wt") as counts_file, open("qc.csv", "wt") as qc_file:
    counts_writer = csv.writer(counts_file)
    qc_writer = csv.writer(qc_file)
    counts_writer.writerow( ['ensembl_id']+ sample_names )
    qc_writer.writerow( ['qc_metric']+ sample_names )
    sorted_genes = sorted(list(sample_dict[sample_names[0]].keys()))
    for gene in sorted_genes:
        output=[gene]
        for sample in sample_names:
            output.append(sample_dict[sample][gene])
        if gene.startswith("N_"):
            qc_writer.writerow(output)
        else:
            counts_writer.writerow(output)

我的slurm脚本：

#!/bin/bash
#SBATCH --job-name=combine_STAR_output
#SBATCH --output=combine_STAR_output.out
#SBATCH --error=combine_STAR_output.err
#SBATCH --nodes=1                         
#SBATCH --ntasks=1  
#SBATCH --cpus-per-task=16
#SBATCH --time=2:00:00
#SBATCH --mem=64G
#SBATCH --partition=ojfat  

python /public/home/linql/STAR-2.7.0e/LNCAP/combine_STAR_output.py

运行结果

Column 3 of LNCAP_Normoxia_S1 stored
Column 3 of LNCAP_Normoxia_S2 stored
Column 3 of PC3_Normoxia_S1 stored
Column 3 of PC3_Normoxia_S2 stored
Column 3 of LNCAP_Hypoxia_S1 stored
Column 3 of LNCAP_Hypoxia_S2 stored
Column 3 of PC3_Hypoxia_S1 stored
Column 3 of PC3_Hypoxia_S2 stored

输出文件的第一列是 Ensembl ID 转录组中每个基因的基因。其余列对应于我们绘制的 8 个样本的每个基因的计数。
将得到以下文件：
在这里插入图片描述

差异表达分析

为了进行差异基因表达分析，我们将使用 R 包。此包提供一组数据专为批量设计的规范化和处理工具 RNA-seq数据和差异基因表达分析。一个很棒的教程可以在此处找到 DESeq2 的细节：https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
在下面的部分中，我将使用 DESeq2 来确定基因 LNCaP 和 PC3 细胞均因缺氧而上调。我还提供自定义函数和可视化功能，简化分析过程并更深入地了解数据。

我们可以使用以下命令从 Bioconductor 安装：DESeq2

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
options(repos=structure(c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")))
BiocManager::install("DESeq2")

library("DESeq2")
library("tidyverse")

现在，我们可以读入 R，并查看数据结构：raw_counts.csv

data <- read.csv("D:\\sra\\raw_counts.csv", header = T, row.names = "ensembl_id")
data <- data[,sort(colnames(data))]
head(data)

运行结果：

                LNCAP_Hypoxia_S1 LNCAP_Hypoxia_S2 LNCAP_Normoxia_S1
ENSG00000000003              608              190               374
ENSG00000000005                0                0                 0
ENSG00000000419             4392             1256              4988
ENSG00000000457              766              204               591
ENSG00000000460              399              131               559
ENSG00000000938                2                2                 2
                LNCAP_Normoxia_S2 PC3_Hypoxia_S1 PC3_Hypoxia_S2
ENSG00000000003               388           1073            333
ENSG00000000005                 0              0              0
ENSG00000000419              5765           3251           1029
ENSG00000000457               653             91             25
ENSG00000000460               509            444            159
ENSG00000000938                 1              0              0
                PC3_Normoxia_S1 PC3_Normoxia_S2
ENSG00000000003             347             964
ENSG00000000005               0               0
ENSG00000000419            1042            2438
ENSG00000000457              31             101
ENSG00000000460             188             471
ENSG00000000938               0               1

我们看到数据有 60676 行和 8 列，其中每行都是一个基因，每列对应于每个基因的reads数每个样本。我们可以通过取列总和：

colSums(data)

LNCAP_Hypoxia_S1  LNCAP_Hypoxia_S2 LNCAP_Normoxia_S1 LNCAP_Normoxia_S2 
        37480174          42634975          40046898          48555048 
 PC3_Hypoxia_S1    PC3_Hypoxia_S2   PC3_Normoxia_S1   PC3_Normoxia_S2 
        37392717          12721451          12754721          37767688

我们观察到每个样本都在不同的“深度”进行了测序，其中每个样本的总读取次数不同。为例如，PC3_Hypoxia_S2只有 1200 万次读取，小于 1/3 的读取次数为 PC3_Hypoxia_S1。这意味着我们不能直接比较样本之间的计数值，无需归一化每列基于其总读取计数。在数据已处理使用，我们将获得可以使用的标准化数据进行比较。DESeq2

创建 DESeq2 对象

要创建 DESeq2 对象，我们可以使用函数，它需要三个项目：DESeqDataSetFromMatrix()

countData- 整数矩阵形式的数据，其中行名对应于基因 ID。这是对象我们从raw_counts.csv
colData- 提供要比较的样本组的表。它还应将每个样本分配给一组（即生物重复）。
design- 一个公式，规定每个基因的计数方式取决于中的变量。我们将使用定义的组在 colData 中。

为了构造每个样本的生物学重复，包括创建基因表达数据矩阵、条件向量和样本信息数据框，以及将这些信息整合到DESeq2数据集中的过程。我们首先定义组和数每组都有样本：colData

# 定义条件：创建一个条件向量，标识每个样本的生物学重复。
condition <- c(rep("LNCAP_Hypoxia", 2), rep("LNCAP_Normoxia", 2), rep("PC3_Hypoxia", 2), rep("PC3_Normoxia", 2))

然后，我们分配各个样本名称（对应于列的名称通过绑定它们到分组类别中在数据帧内：raw_counts.csv

# a创建一个数据框来包含样本信息
my_colData <- as.data.frame(condition)
rownames(my_colData) <- colnames(data)
my_colData

                       condition
LNCAP_Hypoxia_S1   LNCAP_Hypoxia
 LNCAP_Hypoxia_S2   LNCAP_Hypoxia
 LNCAP_Normoxia_S1 LNCAP_Normoxia
LNCAP_Normoxia_S2 LNCAP_Normoxia
 PC3_Hypoxia_S1       PC3_Hypoxia
 PC3_Hypoxia_S2       PC3_Hypoxia
 PC3_Normoxia_S1     PC3_Normoxia
 PC3_Normoxia_S2     PC3_Normoxia

上面命名的向量可用于参数。现在我们已经有了所有必需的输入，我们将创建一个 DESeq2 对象

#创建DESeq2数据集
dds <- DESeqDataSetFromMatrix(countData = data,
                              colData = my_colData,
                              design = ~condition)

然后可以使用函数进行处理，该函数将基于负二项式执行差异表达分析分配。从文档中，此功能将执行：DESeq()

大小因子的估计：estimateSizeFactors
色散估计：estimateDispersions
负二项式 GLM 拟合和 Wald 统计量：nbinomWaldTest

#执行差异表达分析
dds <- DESeq(dds)

estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

dds

class: DESeqDataSet 
dim: 60676 8 
metadata(1): version
assays(4): counts mu H cooks
rownames(60676): ENSG00000000003 ENSG00000000005 ... ENSG00000288587
ENSG00000288588
rowData names(30): baseMean baseVar ... deviance maxCooks
colnames(8): LNCAP_Hypoxia_S1 LNCAP_Hypoxia_S2 ... PC3_Normoxia_S1
PC3_Normoxia_S2
colData names(2): condition sizeFactor

class: 表示该对象的类为 DESeqDataSet，这是 DESeq2 包中的主要数据结构，用于存储差异表达分析的数据。
dim: 表示数据集的维度，包含 60676 个基因（行）和 8 个样本（列）。
metadata: 包含与数据集相关的元数据。这里显示一个元素 version，表示DESeq2包的版本信息。
assays: 包含4个测量值（assays），分别为：counts: 原始计数数据。mu: 估计的期望值（expected values），在拟合模型时使用。H: 计算的散度（dispersion）参数。cooks: Cook’s 距离，用于诊断影响大的样本点。
rownames: 包含基因的名称或标识符，共有 60676 个基因。
rowData names: 包含与基因相关的统计数据的名称，共有 30 项，例如：baseMean: 基础均值，表示所有样本中基因表达水平的平均值。
baseVar: 基础方差，表示基因表达水平的方差。其他如 deviance、maxCooks 等。
colnames: 样本名称，共有 8 个样本。
colData names: 样本相关的变量名称，共有 2 项：
condition: 样本的条件（例如，LNCAP_Hypoxia、LNCAP_Normoxia 等）。
sizeFactor: 样本的大小因子，用于归一化基因表达数据。

我们还可以使用和筛选对象。例如，我们可以通过以下方式找到原始计数数据：dds@$

head(dds@assays@data$counts)

##                 LNCAP_Hypoxia_S1 LNCAP_Hypoxia_S2 LNCAP_Normoxia_S1
## ENSG00000000003              609              693               374
## ENSG00000000005                0                0                 0
## ENSG00000000419             4390             4918              4982
## ENSG00000000457              764              826               593
## ENSG00000000460              398              512               559
## ENSG00000000938                2                2                 2
##                 LNCAP_Normoxia_S2 PC3_Hypoxia_S1 PC3_Hypoxia_S2 PC3_Normoxia_S1
## ENSG00000000003               388           1075            333             348
## ENSG00000000005                 0              0              0               0
## ENSG00000000419              5766           3254           1030            1042
## ENSG00000000457               654             91             25              31
## ENSG00000000460               509            444            160             188
## ENSG00000000938                 1              0              0               0
##                 PC3_Normoxia_S2
## ENSG00000000003             964
## ENSG00000000005               0
## ENSG00000000419            2438
## ENSG00000000457             101
## ENSG00000000460             471
## ENSG00000000938               1

你可以看到，不同样本之间的计数值有很大的差异，这就需要归一化和差异表达分析。

从对象获取原始计数的另一种方法是使用功能。如果您想要归一化计数，我们设置 .例如：

#提取标准化后的计数数据
normalized_counts <- counts(dds, normalized = T)
head(normalized_counts)

##                 LNCAP_Hypoxia_S1 LNCAP_Hypoxia_S2 LNCAP_Normoxia_S1
## ENSG00000000003       332.132920      329.1677896        196.947687
## ENSG00000000005         0.000000        0.0000000          0.000000
## ENSG00000000419      2394.192973     2335.9988300       2623.511696
## ENSG00000000457       416.665930      392.3414058        312.272669
## ENSG00000000460       217.058953      243.1946728        294.368334
## ENSG00000000938         1.090749        0.9499792          1.053196
##                 LNCAP_Normoxia_S2 PC3_Hypoxia_S1 PC3_Hypoxia_S2 PC3_Normoxia_S1
## ENSG00000000003       165.0244594      1375.2607     1186.97279      1061.85621
## ENSG00000000005         0.0000000         0.0000        0.00000         0.00000
## ENSG00000000419      2452.3995699      4162.8821     3671.41735      3179.46601
## ENSG00000000457       278.1597847       116.4174       89.11207        94.59064
## ENSG00000000460       216.4882728       568.0146      570.31726       573.64646
## ENSG00000000938         0.4253208         0.0000        0.00000         0.00000
##                 PC3_Normoxia_S2
## ENSG00000000003     1113.540719
## ENSG00000000005        0.000000
## ENSG00000000419     2816.195304
## ENSG00000000457      116.667648
## ENSG00000000460      544.063982
## ENSG00000000938        1.155125

如您所见，归一化计数的值与正如预期的那样，原始计数。不同样本间的总读取数差异被消除，这些标准化计数可用于比较样本之间的基因表达值。

从 BioMart 制作注释文件

目前，这些基因使用其 Ensembl ID 进行标记（例如， ENSG00000111640）。如果能翻译每个 Ensembl ID，那就太好了添加到其基因名称中，这样我们就可以搜索到我们的基因的数据兴趣，而无需每次都查找 Ensembl ID。我们会的在生成时，还需要知道每个 Ensembl ID 的基因名称可视化和差异基因列表。
为了将 Ensembl ID 映射到它们的基因名称，我们必须构建一个在 ensembl.org 使用 BioMart 的基因组注释文件。操作方法如下这样做：

前往 https://www.ensembl.org/biomart/martview/c59a1d3a9024005e61d1d93feb14ddf7 的 Ensembl BioMart
选择要使用的数据库，在本例中为“Ensembl Genes 99”。
选择要使用的注释，在本例中为“人类基因” （GRCh38.p14）”。
在左侧边栏上，单击“属性”，这将允许您选择要映射的项目。对于此分析，您需要：
“基因稳定 ID”，即我们raw_counts.csv中的 Ensembl ID 文件（例如ENSG00000111640）
“基因名”，这是最常用的基因名称（例如：GAPDH）
“基因类型”，指定哪些基因是蛋白质编码的。如果我们想进行基因集富集分析（GSEA），我们会需要这个类别。
您可以取消勾选“基因稳定 ID 版本”、“转录本稳定” ID“和”成绩单稳定 ID 版本”。
单击结果按钮，选中“仅唯一结果”，选择 “导出为 csv”，然后点击“开始”按钮！
批注应以“mart_export.txt”的形式下载到您的计算机上。我们可以将其重命名为“GRCh38.p13_annotation.csv”，并将其读入 R：

annotation <- read.csv("D:\\sra\\GRCh38.p14_annotation.csv", header = T, stringsAsFactors = F)
head(annotation)

    Gene.stable.ID Gene.name      Gene.type
1 ENSG00000210049     MT-TF        Mt_tRNA
2 ENSG00000211459   MT-RNR1        Mt_rRNA
3 ENSG00000210077     MT-TV        Mt_tRNA
4 ENSG00000210082   MT-RNR2        Mt_rRNA
5 ENSG00000209082    MT-TL1        Mt_tRNA
6 ENSG00000198888    MT-ND1 protein_coding

我们观察到有三列对应于注释我们选择的。要将基因名称映射到另一个数据集，我们可以联接两个基于 Ensembl ID 的表，使用 Ensembl ID 中的联接函数。例如，我可以使用right_join()归一化将注释映射到Counts 表：

install.packages("tibble")
library(tibble)
install.packages("dplyr")
library(dplyr)
normalized_counts <- rownames_to_column(as.data.frame(normalized_counts), var = "ensembl_id")
annotated_data <- right_join(annotation, normalized_counts, by = c("Gene.stable.ID" = "ensembl_id"))
head(annotated_data)

   Gene.stable.ID Gene.name      Gene.type LNCAP_Hypoxia_S1 LNCAP_Hypoxia_S2
1 ENSG00000210049     MT-TF        Mt_tRNA     1.090749e+00     2.849938e+00
2 ENSG00000211459   MT-RNR1        Mt_rRNA     1.970219e+04     2.040840e+04
3 ENSG00000210077     MT-TV        Mt_tRNA     0.000000e+00     0.000000e+00
4 ENSG00000210082   MT-RNR2        Mt_rRNA     1.560468e+05     1.578694e+05
5 ENSG00000209082    MT-TL1        Mt_tRNA     3.817620e+00     4.749896e+00
6 ENSG00000198888    MT-ND1 protein_coding     5.395278e+04     5.238375e+04
  LNCAP_Normoxia_S1 LNCAP_Normoxia_S2 PC3_Hypoxia_S1 PC3_Hypoxia_S2
1      2.422351e+01          11.90898   1.407243e+02   1.461438e+02
2      4.357020e+04       46733.82108   6.888532e+06   6.490617e+06
3      2.106392e+00           0.00000   1.791037e+02   1.461438e+02
4      5.286887e+05      443639.33738   1.979179e+07   1.883179e+07
5      1.685114e+01          25.09393   4.989318e+01   8.911207e+01
6      1.509241e+05      178563.27024   4.158532e+04   3.675695e+04
 PC3_Normoxia_S1 PC3_Normoxia_S2
1    5.797491e+01    1.594073e+02
2    6.654287e+06    7.287938e+06
3    2.410536e+02    3.881221e+02
4    1.456696e+07    1.793767e+07
5    5.187229e+01    1.189779e+02
6    6.430638e+04    8.450896e+04

write.csv(annotated_data, file = "D:\\sra\\gene_annotated_normalized_counts.csv")

现在，基因名称直接放在 Ensembl ID 旁边。我们可以导出此规范化计数文件并将其提供给协作人员科学家们，他们将能够快速搜索他们的基因利息！有关如何绘制给定的归一化计数的示例跨样本的基因，请参阅以下部分：DE 的可视化结果。

样品变异性的可视化

此时，我们可以开始评估样本间的变异性。这对于确定数据集中是否有异常值非常重要，以及您的数据是否有意义。例如，基因生物/技术重复的表达模式应观察彼此非常相似。在本节中，我们将介绍一些内容可用于此目的的可视化效果。为了使用这些函数中，我们首先对数据使用：vst()

vsd <- vst(dds, blind = TRUE)

vst 函数用于对 DESeqDataSet 对象进行方差稳定化变换（variance stabilizing transformation），这在差异表达分析和下游分析（如聚类和可视化）中非常有用。blind = TRUE: 表示在执行变换时忽略实验设计。这通常用于探索性数据分析。

距离图
距离图计算每个单独样本的表达式值。这是一种识别方法哪些样本在基因表达方面最密切相关。

install.packages("pheatmap")
library(pheatmap)

plotDists = function (vsd.obj) {
  sampleDists <- dist(t(assay(vsd.obj)))# 计算样本之间的距离
  sampleDistMatrix <- as.matrix( sampleDists )# 将距离对象转换为矩阵
  rownames(sampleDistMatrix) <- paste( vsd.obj$condition )# 设置矩阵的行名为样本的条件
  colors <- colorRampPalette( rev(RColorBrewer::brewer.pal(9, "Blues")) )(255)
  pheatmap::pheatmap(sampleDistMatrix,
           clustering_distance_rows = sampleDists,
           clustering_distance_cols = sampleDists,
           col = colors)
}
plotDists(vsd)

在这里插入图片描述

我们观察到所有 LNCaP 样本聚集在一起，所有 PC3 样本聚集在一起。在每种细胞类型中，缺氧样本聚集在一起，常氧样品聚集在一起。这些聚类模式是有道理的，表明实验是成功的。

可变基因热图
我们可以确定两者之间聚类的基因采样并将它们显示在热图中。为此，我们首先获得顶部通过选择具有最高值。然后，我们插入这些的归一化计数使用将基因转换为热图。热图将根据表达相似度和显示对样本进行聚类与右侧每个簇相关的基因。请注意，这些基因只是根据跨领域的变异性来选择的样品。不执行任何统计检验来确定它们是否实际上，两组之间有显著差异。rowVars()pheatmap::pheatmap()

variable_gene_heatmap <- function (vsd.obj, num_genes = 500, annotation, title = "") {
  brewer_palette <- "RdBu"
  # 创建颜色梯度
  ramp <- colorRampPalette( RColorBrewer::brewer.pal(11, brewer_palette))
  mr <- ramp(256)[256:1]
  # 获取标准化后的计数数据
  stabilized_counts <- assay(vsd.obj)
  # 计算每行（基因）的方差
  row_variances <- rowVars(stabilized_counts)
  # 获取方差最大的前 num_genes 个基因
  top_variable_genes <- stabilized_counts[order(row_variances, decreasing=T)[1:num_genes],]
  # 减去每行的均值，使得每个基因的均值为零
  top_variable_genes <- top_variable_genes - rowMeans(top_variable_genes, na.rm=T)
  # 使用基因名称替换 Ensembl ID
  gene_names <- annotation$Gene.name[match(rownames(top_variable_genes), annotation$Gene.stable.ID)]
  rownames(top_variable_genes) <- gene_names
  # 构建没有 sizeFactors 的 colData 以用于热图标签
  coldata <- as.data.frame(vsd.obj@colData)
  coldata$sizeFactor <- NULL
  pheatmap::pheatmap(top_variable_genes, color = mr, annotation_col = coldata, fontsize_col = 8, fontsize_row = 250/num_genes, border_color = NA, main = title)
}

variable_gene_heatmap(vsd, num_genes = 40, annotation = annotation)

在这里插入图片描述
我们看到聚类与我们在距离图。有一组独特的基因是特异性的 LNCaP 或 PC3 细胞。鉴于我们正在分析来自单元格的数据线，聚类定义如此明确也就不足为奇了。如果您正在处理高度异质的样品（例如组织），样本之间可能存在显着差异。了解基因定义样本之间的差异对于评估是否存在样品制备问题，例如细胞质量差纯度或组织污染。

PCA图
主成分图是观察样本之间的总体相似性和差异性。与往常一样，每个最相似的样本其他的就其基因表达值而言将更接近每个其他在情节上。此图使用笛卡尔坐标系，并且显示的轴对应于前两个主分量这解释了数据中的大部分可变性。

install.packages("ggplot2")
library(ggplot2)

plot_PCA = function (vsd.obj) {
  pcaData <- plotPCA(vsd.obj,  intgroup = c("condition"), returnData = T)
  percentVar <- round(100 * attr(pcaData, "percentVar"))
  ggplot(pcaData, aes(PC1, PC2, color=condition)) +
    geom_point(size=3) +
    labs(x = paste0("PC1: ",percentVar[1],"% variance"),
         y = paste0("PC2: ",percentVar[2],"% variance"),
         title = "PCA Plot colored by condition") +
    ggrepel::geom_text_repel(aes(label = name), color = "black")
}
plot_PCA(vsd)

在这里插入图片描述

我们观察到 LNCaP 样本与 PC3 相距最远 x 轴上的样本，这解释了高达 99% 的数据集中的方差。在每种细胞类型中，常氧样品是从缺氧样品中分离出来。引入的可变性量与缺氧相比，缺氧的固有差异较小 LNCaP 和 PC3 细胞系。如果我们对缺氧的影响感兴趣在基因表达方面，我们应该将 LNCaP 样本与 PC3 分开。首先进行采样，创建单独的 DESeq2 对象，然后确定缺氧对每个细胞系的影响单独存在。我将介绍如何在下一节中执行此操作。

创建子集 DESeq2 对象

因为数据中的很多可变性是由于细胞系的背景，我们必须创建单独的 DESeq2 对象用于 LNcaP 和 PC3 样品。下面的函数将使用您的原始计数数据和两者创建 DESeq2 对象您指定进行比较的组。它将从中选择列与组名称匹配的原始数据，根据您的组名称构建矩阵，并生成使用和的 DESeq2 对象。

generate_DESeq_object <- function (my_data, groups) {
  data_subset1 <- my_data[,grep(str_c("^", groups[1]), colnames(my_data))]
  data_subset2 <- my_data[,grep(str_c("^", groups[2]), colnames(my_data))]
  my_countData <- cbind(data_subset1, data_subset2)
  condition <- c(rep(groups[1],ncol(data_subset1)), rep(groups[2],ncol(data_subset2)))
  my_colData <- as.data.frame(condition)
  rownames(my_colData) <- colnames(my_countData)
  print(my_colData)
  dds <- DESeqDataSetFromMatrix(countData = my_countData,
                              colData = my_colData,
                              design = ~ condition)
  dds <- DESeq(dds, quiet = T)
  return(dds)
}

该函数应打印用于生成每个 DESeq2 对象，因此我们可以仔细检查以确保适当的已选择样本列，并且分组正确。colData

lncap <- generate_DESeq_object(data, c("LNCAP_Hypoxia", "LNCAP_Normoxia"))

                       condition
LNCAP_Hypoxia_S1   LNCAP_Hypoxia
LNCAP_Hypoxia_S2   LNCAP_Hypoxia
LNCAP_Normoxia_S1 LNCAP_Normoxia
LNCAP_Normoxia_S2 LNCAP_Normoxia

pc3 <- generate_DESeq_object(data, c("PC3_Hypoxia", "PC3_Normoxia"))

                   condition
PC3_Hypoxia_S1   PC3_Hypoxia
PC3_Hypoxia_S2   PC3_Hypoxia
PC3_Normoxia_S1 PC3_Normoxia
PC3_Normoxia_S2 PC3_Normoxia

现在我们已经分离了 LNCaP 和 PC3 样本，应该显示富集缺氧与常氧条件。

lncap_vsd <- vst(lncap, blind = T)
pc3_vsd <- vst(pc3, blind = T)
a <- variable_gene_heatmap(lncap_vsd, 30, annotation = annotation, title = "LNCaP variable genes")

b <- variable_gene_heatmap(pc3_vsd, 30, annotation = annotation, title = "PC3 variable genes")

gridExtra::grid.arrange(a[[4]],b[[4]], nrow = 1)

在这里插入图片描述

这是一组与我们观察到的基因截然不同的基因将两种细胞系的样品合并。我们的差异基因表达分析现在应该产生更深入的基因我们关心的比较：缺氧与常氧。

提取 DE 结果

现在，我们已经评估了样本间的变异性并获得了确信我们的实验是正确执行的，我们已经准备好了开始回答最紧迫的问题：什么是基因每个细胞中因缺氧而上调或下调差异线？在本节中，我将解释如何生成微分来自 DESeq2 对象的表达基因列表。

为了从 DESeq2 对象中提取差异表达的基因，我们将使用以下函数：results()

results(lncap, contrast = c("condition", "LNCAP_Hypoxia", "LNCAP_Normoxia"))

log2 fold change (MLE): condition LNCAP_Hypoxia vs LNCAP_Normoxia 
Wald test p-value: condition LNCAP_Hypoxia vs LNCAP_Normoxia 
DataFrame with 60676 rows and 6 columns
                 baseMean log2FoldChange     lfcSE      stat
                <numeric>      <numeric> <numeric> <numeric>
ENSG00000000003   376.956       0.911073  0.196795  4.629561
ENSG00000000005     0.000             NA        NA        NA
ENSG00000000419  3536.567      -0.108109  0.112746 -0.958868
ENSG00000000457   499.550       0.439211  0.177746  2.471009
ENSG00000000460   350.256      -0.160835  0.217695 -0.738810
...                   ...            ...       ...       ...
ENSG00000288584    0.0000             NA        NA        NA
ENSG00000288585   52.5050      -1.016896  0.480612  -2.11584
ENSG00000288586   55.8962      -0.703165  0.451981  -1.55574
ENSG00000288587    0.0000             NA        NA        NA
ENSG00000288588    0.0000             NA        NA        NA
                     pvalue        padj
                  <numeric>   <numeric>
ENSG00000000003 3.66441e-06 3.83609e-05
ENSG00000000005          NA          NA
ENSG00000000419 3.37625e-01 5.12463e-01
ENSG00000000457 1.34733e-02 4.44083e-02
ENSG00000000460 4.60022e-01 6.28210e-01
...                     ...         ...
ENSG00000288584          NA          NA
ENSG00000288585   0.0343587   0.0938579
ENSG00000288586   0.1197702   0.2432893
ENSG00000288587          NA          NA
ENSG00000288588          NA          NA

我们观察到，将缺氧条件与常氧条件进行了比较条件，并且存在与统计信息相对应的列为每个基因计算。您应该支付最多的列注意的是、和。

所有样本中基因表达的平均值。我们可以评估它是否是一个高表达的基因，或者它是否只是一个勉强高于基线。baseMean
这些指标衡量了基因在缺氧条件与常氧条件相比。阴性值表示常氧组的表达量较高。log2FoldChange
log2FoldChange 的标准误。lfcSE
Wald 统计量。 stat
Wald 检验的 p 值。pvalue
多重检验校正后的 p 值（例如，Benjamini-Hochberg 校正）。该列是衡量基因可能性的指标被比较的两种条件之间存在显着差异。padj

我们将使用这些指标来过滤、组织和导出以下列表差异表达基因。我选择使用以下条件：

根据、选择显著不同的基因，在低于指定的截止值。padjpadj
排序方式以识别出含量最高的基因向上或向下。log2FoldChange
（可选）过滤掉表达非常低（低值）的基因，这些基因在生物学上不太可能是相关。为此，我们可以使用或计算另一个指标称为“百万计数”（CPM）。每个计数百万值更直观，因为我们可以使用类似的截止值在不同测序实验中的价值。baseMeanbaseMean
生成一个排名基因列表文件（），它是所有基因列表的列表数据集中的基因按排序。此文件是基因集富集分析所必需的，我将对此进行解释稍后在 GSEA 部分。.rnklog2FoldChange

由于对数据进行排序/过滤有很多种变化，所以我写了执行多个过滤步骤的函数并以 CSV 文件的形式提供输出。我导出这些文件来提取和注释差异表达基因，以便它们可以很容易地被其他实验室科学家分享和解释而无需将数据读入 R 中。generate_DE_results()

generate_DE_results <- function (dds, comparisons, padjcutoff = 0.001, log2cutoff = 0.5, cpmcutoff = 2) {
  # 从原始计数数据计算每百万次计数（CPM）的平均值。
  raw_counts <- counts(dds, normalized = F)
  cpms <- enframe(rowMeans(edgeR::cpm(raw_counts)))
  colnames(cpms) <- c("ensembl_id", "avg_cpm")
  
  # 提取指定对比的DESeq结果，删除不必要的列。
  res <- results(dds, contrast = c("condition", comparisons[1], comparisons[2]))[,-c(3,4)]
  
  # 将基因注释和CPM数据与DESeq结果结合。
  res <- as_tibble(res, rownames = "ensembl_id")
  # read in the annotation and append it to the data
  my_annotation <- read.csv("GRCh38.p14_annotation.csv", header = T, stringsAsFactors = F)
  res <- left_join(res, my_annotation, by = c("ensembl_id" = "Gene.stable.ID"))
  # append the average cpm value to the results data
  res <- left_join(res, cpms, by = c("ensembl_id" = "ensembl_id"))
  
  # 合并标准化的计数数据，并按log2FoldChange排序。
  normalized_counts <- round(counts(dds, normalized = TRUE),3)
  pattern <- str_c(comparisons[1], "|", comparisons[2])
  combined_data <- as_tibble(cbind(res, normalized_counts[,grep(pattern, colnames(normalized_counts))] ))
  combined_data <- combined_data[order(combined_data$log2FoldChange, decreasing = T),]
  
  # 生成GSEA排名文件, selecting only protein coding genes
  res_prot <- res[which(res$Gene.type == "protein_coding"),]
  res_prot_ranked <- res_prot[order(res_prot$log2FoldChange, decreasing = T),c("Gene.name", "log2FoldChange")]
  res_prot_ranked <- na.omit(res_prot_ranked)
  res_prot_ranked$Gene.name <- str_to_upper(res_prot_ranked$Gene.name)
  
  #生成排序列表 with the indicated cutoff values
  res <- res[order(res$log2FoldChange, decreasing=TRUE ),]
  de_genes_padj <- res[which(res$padj < padjcutoff),]
  de_genes_log2f <- res[which(abs(res$log2FoldChange) > log2cutoff & res$padj < padjcutoff),]
  de_genes_cpm <- res[which(res$avg_cpm > cpmcutoff & res$padj < padjcutoff),]
  
  # 输出
  write.csv (de_genes_padj, file = paste0(comparisons[1], "_vs_", comparisons[2], "_padj_cutoff.csv"), row.names =F)
  write.csv (de_genes_log2f, file = paste0(comparisons[1], "_vs_", comparisons[2], "_log2f_cutoff.csv"), row.names =F)
  write.csv (de_genes_cpm, file = paste0(comparisons[1], "_vs_", comparisons[2], "_cpm_cutoff.csv"), row.names =F)
  write.csv (combined_data, file = paste0(comparisons[1], "_vs_", comparisons[2], "_allgenes.csv"), row.names =F)
  write.table (res_prot_ranked, file = paste0(comparisons[1], "_vs_", comparisons[2], "_rank.rnk"), sep = "\t", row.names = F, quote = F)
  
  writeLines( paste0("For the comparison: ", comparisons[1], "_vs_", comparisons[2], ", out of ", nrow(combined_data), " genes, there were: \n", 
               nrow(de_genes_padj), " genes below padj ", padjcutoff, "\n",
               nrow(de_genes_log2f), " genes below padj ", padjcutoff, " and above a log2FoldChange of ", log2cutoff, "\n",
               nrow(de_genes_cpm), " genes below padj ", padjcutoff, " and above an avg cpm of ", cpmcutoff, "\n",
               "Gene lists ordered by log2fchange with the cutoffs above have been generated.") )
  gene_count <- tibble (cutoff_parameter = c("padj", "log2fc", "avg_cpm" ), 
                        cutoff_value = c(padjcutoff, log2cutoff, cpmcutoff), 
                        signif_genes = c(nrow(de_genes_padj), nrow(de_genes_log2f), nrow(de_genes_cpm)))
  invisible(gene_count)
}

运行该函数后，控制台应显示通过每个过滤的基因：

lncap_output <- generate_DE_results (lncap, c("LNCAP_Hypoxia", "LNCAP_Normoxia"))

For the comparison: LNCAP_Hypoxia_vs_LNCAP_Normoxia, out of 60676 genes, there were: 
2448 genes below padj 0.001
2298 genes below padj 0.001 and above a log2FoldChange of 0.5
2350 genes below padj 0.001 and above an avg cpm of 2
Gene lists ordered by log2fchange with the cutoffs above have been generated.

pc3_output <- generate_DE_results(pc3, c("PC3_Hypoxia", "PC3_Normoxia"))

For the comparison: PC3_Hypoxia_vs_PC3_Normoxia, out of 60676 genes, there were: 
1608 genes below padj 0.001
1359 genes below padj 0.001 and above a log2FoldChange of 0.5
1564 genes below padj 0.001 and above an avg cpm of 2
Gene lists ordered by log2fchange with the cutoffs above have been generated.

我们还应该看到一组 csv 文件显示在与 R 脚本。我们可以读取输出文件并观察文件结构是什么样的：read.csv()

res <- read.csv("LNCAP_Hypoxia_vs_LNCAP_Normoxia_allgenes.csv", header = T)
head(res)

       ensembl_id  baseMean log2FoldChange       pvalue         padj
1 ENSG00000283339  28.69060       8.824554 1.437852e-08 2.432106e-07
2 ENSG00000203664  26.76639       8.746350 2.306888e-08 3.799485e-07
3 ENSG00000261051  13.74824       7.776062 3.195089e-06 3.401621e-05
4 ENSG00000276107  13.03045       7.640454 7.254747e-06 7.012870e-05
5 ENSG00000171954  11.49608       7.550623 1.151999e-05           NA
6 ENSG00000105329 146.91475       7.422230 3.254435e-26 2.761015e-24
                           Gene.type Gene.name   avg_cpm
1             unprocessed_pseudogene           0.9822830
2 transcribed_unprocessed_pseudogene    OR2W5P 0.9141435
3                             lncRNA           0.4700321
4                             lncRNA THBS1-IT1 0.4484480
5                     protein_coding   CYP4F22 0.3914903
6                     protein_coding     TGFB1 5.0196594
  LNCAP_Hypoxia_S1 LNCAP_Hypoxia_S2 LNCAP_Normoxia_S1 LNCAP_Normoxia_S2
1           53.855           60.907             0.000             0.000
2           59.399           47.667             0.000             0.000
3           28.511           26.482             0.000             0.000
4           15.048           37.074             0.000             0.000
5           30.095           15.889             0.000             0.000
6          311.250          272.760             3.032             0.618

该文件将在下一节用于可视化。它是我们精心策划的所有数据的汇总 FAR：输出、基因注释和归一化每个基因都很重要。我们可以进一步组织/过滤数据基于用户偏好的可视化。其他 csv 文件（例如）可以分发给仅对基因列表感兴趣。

DE 结果的可视化效果

现在我们已经从 DESeq2 对象中提取了结果，我们可以开始可视化差异表达的基因。在本节中，我将展示我为回答问题而编写的一些绘图函数如：

给定的兴趣基因如何在各组之间表达？
顶部微分的样本间变异性是多少表达基因？
是否存在更多上调或下调的基因？
如果我们进行多次比较，这是不同的表达的基因是共享的？
哪些基因通路在一个群体中富集而在另一个群体中富集？

PlotCounts（升级）
我们通常对绘制某些归一化计数感兴趣感兴趣的基因，以便了解数据的传播和基因在每个组中的表达方式。这尤其是对在我们的 DESeq2 结果。绘制每个样本的归一化计数的一种方法是通过使用内置函数。然而，这功能还有很多不足之处。举个例子，我绘制了基因 IGFBP1（一种缺氧诱导基因）使用其下面的 Ensembl ID。DESeq2::plotCounts()

plotCounts(dds, gene="ENSG00000146678", intgroup="condition")

在这里插入图片描述

如您所见，如果样本名称及其名字太长。我们可以改用 ggplot2 来绘制数据，通过以下方式使用以下方法提取归一化计数：returnData = TRUE

d <- plotCounts(dds, gene="ENSG00000146678", intgroup="condition", returnData=TRUE)
d

                      count      condition
LNCAP_Hypoxia_S1   0.500000  LNCAP_Hypoxia
LNCAP_Hypoxia_S2   0.500000  LNCAP_Hypoxia
LNCAP_Normoxia_S1  0.500000 LNCAP_Normoxia
LNCAP_Normoxia_S2  0.500000 LNCAP_Normoxia
PC3_Hypoxia_S1    25.793111    PC3_Hypoxia
PC3_Hypoxia_S2    22.031789    PC3_Hypoxia
PC3_Normoxia_S1    0.500000   PC3_Normoxia
PC3_Normoxia_S2    4.448567   PC3_Normoxia

ggplot(d, aes(x=condition, y=count)) + 
  geom_point(position=position_jitter(w=0.1,h=0), aes (color = condition)) + 
  scale_y_log10(breaks=c(25,100,400))

在这里插入图片描述
然而，这仍然非常令人沮丧，因为我们必须这样做输入 Ensembl ID 而不是基因名称。研究人员通常希望按基因名称搜索，此功能迫使我们查找每次都接收 ID。为了解决这个问题，我编写了自己的函数来接受基因名称、Ensembl ID 或索引如。我还提供了两种方法来规范化数据：按百万分之一（CPM）或使用内置的 DESeq2 标准化计数。plot_counts()which.min(res$padj)

plot_counts <- function (dds, gene, normalization = "DESeq2"){
  annotation <- read.csv("GRCh38.p14_annotation.csv", header = T, stringsAsFactors = F)
  # 获取标准化数据
  if (normalization == "cpm") {
    normalized_data <- cpm(counts(dds, normalized = F)) # 按每百万次计数标准化（counts per million）
  } else if (normalization == "DESeq2")
    normalized_data <- counts(dds, normalized = T) # 使用DESeq2 计数标准化
  #  获取样本组信息
  condition <- dds@colData$condition
  # 确定基因ID
  if (is.numeric(gene)) { # 检查是否提供了索引或ENSEMBL ID
    if (gene%%1==0 )
      ensembl_id <- rownames(normalized_data)[gene]
    else
      stop("Invalid index supplied.")
  } else if (gene %in% annotation$Gene.name){ # 检查是否提供了基因名
    ensembl_id <- annotation$Gene.stable.ID[which(annotation$Gene.name == gene)]
  } else if (gene %in% annotation$Gene.stable.ID){
    ensembl_id <- gene
  } else {
    stop("Gene not found. Check spelling.")
  }
  # 获取基因表达数据
  expression <- normalized_data[ensembl_id,]
  gene_name <- annotation$Gene.name[which(annotation$Gene.stable.ID == ensembl_id)]
  # 构建数据框并绘制图形
  gene_tib <- tibble(condition = condition, expression = expression)
  ggplot(gene_tib, aes(x = condition, y = expression))+
    geom_boxplot(outlier.size = NULL)+
    geom_point()+
    labs (title = paste0("Expression of ", gene_name, " - ", ensembl_id), x = "group", y = paste0("Normalized expression (", normalization , ")"))+
    theme(axis.text.x = element_text(size = 11), axis.text.y = element_text(size = 11))
}

plot_counts(dds, "IGFBP1")

在这里插入图片描述

我们观察到 IGFBP1 在缺氧条件下升高相对于相应的常氧对照。有趣的是， PC3 细胞中 IGFBP1 的上调远高于 LNCaP 细胞。这是一个很好的例子，说明这个图如何用于确定相对样本之间的表达水平。

差异基因热图谱
用于显示 RNA-seq 结果的最常见热图是差异基因热图。这是一个热图，其中归一化值对于每个样本，绘制了上调最深的基因，这些基因是对照组间差异显著。从 results 对象，我们首先过滤结果，然后取具有最大价值的基因。归一化每个基因的计数按行（跨样本）缩放，以便强调对照组之间的差异（缺氧 vs 常氧）。此可视化效果的目的是让您了解如何实现可变数据从一个样本到另一个样本，以及快速检查排名靠前的差异基因。我编写了函数来生成此图，并与参数一起使用。

res <- read.csv ("LNCAP_Hypoxia_vs_LNCAP_Normoxia_allgenes.csv", header = T)

DE_gene_heatmap <- function(res, padj_cutoff = 0.0001, ngenes = 20) {
  # generate the color palette
  brewer_palette <- "RdBu"
  ramp <- colorRampPalette(RColorBrewer::brewer.pal(11, brewer_palette))
  mr <- ramp(256)[256:1]
  # obtain the significant genes and order by log2FoldChange
  significant_genes <- res %>% filter(padj < padj_cutoff) %>% arrange (desc(log2FoldChange)) %>% head (ngenes)
  heatmap_values <- as.matrix(significant_genes[,-c(1:8)])
  rownames(heatmap_values) <- significant_genes$Gene.name
  # plot the heatmap using pheatmap
  pheatmap::pheatmap(heatmap_values, color = mr, scale = "row", fontsize_col = 10, fontsize_row = 200/ngenes, fontsize = 5, border_color = NA)
}
DE_gene_heatmap(res, 0.001, 30)

以下是作者原图
在这里插入图片描述
以下是我生成的

我们观察到差异基因的表达非常一致在组内。我们还能够确定LNCAP_Hypoxia_S1具有前~10个基因的表达高于LNCAP_Hypoxia_S2。

火山图

批量 RNA-seq 的另一个常见可视化是火山图，它是一个散点图，显示比较中每个基因的值与值的关系。我们可以使用这个图为了快速识别高度上调或下调的基因，或识别具有高度显着 P 值的基因。下面是一个显示缺氧中最重要基因的火山图示例与LNCaP细胞的常氧比较。log2FoldChange-log(padj)

res <- read.csv ("LNCAP_Hypoxia_vs_LNCAP_Normoxia_allgenes.csv", header = T)

plot_volcano <- function (res, padj_cutoff, nlabel = 10, label.by = "padj"){
  # assign significance to results based on padj
  res <- mutate(res, significance=ifelse(res$padj<padj_cutoff, paste0("padj < ", padj_cutoff), paste0("padj > ", padj_cutoff)))
  res = res[!is.na(res$significance),]
  significant_genes <- res %>% filter(significance == paste0("padj < ", padj_cutoff))
  
  # get labels for the highest or lowest genes according to either padj or log2FoldChange
  if (label.by == "padj") {
    top_genes <- significant_genes %>% arrange(padj) %>% head(nlabel)
    bottom_genes <- significant_genes %>% filter (log2FoldChange < 0) %>% arrange(padj) %>% head (nlabel)
  } else if (label.by == "log2FoldChange") {
    top_genes <- head(arrange(significant_genes, desc(log2FoldChange)),nlabel)
    bottom_genes <- head(arrange(significant_genes, log2FoldChange),nlabel)
  } else
    stop ("Invalid label.by argument. Choose either padj or log2FoldChange.")
  
  ggplot(res, aes(log2FoldChange, -log(padj))) +
    geom_point(aes(col=significance)) + 
    scale_color_manual(values=c("red", "black")) + 
    ggrepel::geom_text_repel(data=top_genes, aes(label=head(Gene.name,nlabel)), size = 3)+
    ggrepel::geom_text_repel(data=bottom_genes, aes(label=head(Gene.name,nlabel)), color = "#619CFF", size = 3)+
    labs ( x = "Log2FoldChange", y = "-(Log normalized p-value)")+
    geom_vline(xintercept = 0, linetype = "dotted")+
    theme_minimal()
}

plot_volcano(res, 0.0005, nlabel = 15, label.by = "padj")

在这里插入图片描述

与下调的基因相比，上调的基因被标记为不同的颜色。我们观察到，还有更多更重要的东西上调的基因比下调的基因多。鉴于我们是寻找在缺氧条件下生长的细胞中的基因表达变化与常氧条件相比，这表明细胞可能正在急剧增加反应通路的向上转录。

4.LogFoldChange 比较图
鉴于我们已经进行了差异基因表达分析两种不同的细胞系，我们想比较这两种细胞系它们的基因表达如何对缺氧做出反应。我们想要回答以下问题：

缺氧诱导的基因在两个细胞中是否相似线？
两种细胞系共享哪些重要基因，以及哪些每个细胞系的基因都是独一无二的吗？
为此，我们可以从缺氧与常氧中提取重要基因比较两条线，并使用下面的自定义函数生成其值的散点图。此功能将显示该基因按颜色列出编码两者共享的重要基因细胞系，以及每个细胞独特的细胞系线。

res1 <- read.csv ("LNCAP_Hypoxia_vs_LNCAP_Normoxia_allgenes.csv", header = T)
res2 <- read.csv ("PC3_Hypoxia_vs_PC3_Normoxia_allgenes.csv", header = T)

compare_significant_genes <- function (res1, res2, padj_cutoff=0.0001, ngenes=250, nlabel=10, samplenames=c("comparison1", "comparison2"), title = "" ) {
  # get list of most upregulated or downregulated genes for each results table
  genes1 <- rbind(head(res1[which(res1$padj < padj_cutoff),], ngenes), tail(res1[which(res1$padj < padj_cutoff),], ngenes))
  genes2 <- rbind(head(res2[which(res2$padj < padj_cutoff),], ngenes), tail(res2[which(res2$padj < padj_cutoff),], ngenes))
 
   # combine the data from both tables
  de_union <- union(genes1$ensembl_id,genes2$ensembl_id)
  res1_union <- res1[match(de_union, res1$ensembl_id),][c("ensembl_id", "log2FoldChange", "Gene.name")]
  res2_union <- res2[match(de_union, res2$ensembl_id),][c("ensembl_id", "log2FoldChange", "Gene.name")]
  combined <- left_join(res1_union, res2_union, by = "ensembl_id", suffix = samplenames )
  
  # identify overlap between genes in both tables
  combined$de_condition <- 1 # makes a placeholder column
  combined$de_condition[which(combined$ensembl_id %in% intersect(genes1$ensembl_id,genes2$ensembl_id))] <- "Significant in Both"
  combined$de_condition[which(combined$ensembl_id %in% setdiff(genes1$ensembl_id,genes2$ensembl_id))] <- paste0("Significant in ", samplenames[1])
  combined$de_condition[which(combined$ensembl_id %in% setdiff(genes2$ensembl_id,genes1$ensembl_id))] <- paste0("Significant in ", samplenames[2])
  combined[is.na(combined)] <- 0
  
  # find the top most genes within each condition to label on the graph
  label1 <- rbind(head(combined[which(combined$de_condition==paste0("Significant in ", samplenames[1])),],nlabel),
                  tail(combined[which(combined$de_condition==paste0("Significant in ", samplenames[1])),],nlabel))
  label2 <- rbind(head(combined[which(combined$de_condition==paste0("Significant in ", samplenames[2])),],nlabel),
                  tail(combined[which(combined$de_condition==paste0("Significant in ", samplenames[2])),],nlabel))
  label3 <- rbind(head(combined[which(combined$de_condition=="Significant in Both"),],nlabel),
                  tail(combined[which(combined$de_condition=="Significant in Both"),],nlabel))
  combined_labels <- rbind(label1,label2,label3)
  
  # plot the genes based on log2FoldChange, color coded by significance
  ggplot(combined, aes_string(x = paste0("log2FoldChange", samplenames[1]), y = paste0("log2FoldChange", samplenames[2]) )) +
      geom_point(aes(color = de_condition), size = 0.7)+
      scale_color_manual(values= c("#00BA38", "#619CFF", "#F8766D"))+
      ggrepel::geom_text_repel(data= combined_labels, aes_string(label=paste0("Gene.name", samplenames[1]), color = "de_condition"), show.legend = F, size=3)+
      geom_vline(xintercept = c(0,0), size = 0.3, linetype = 2)+ 
      geom_hline(yintercept = c(0,0), size = 0.3, linetype = 2)+
      labs(title = title,x = paste0("log2FoldChange in ", samplenames[1]), y = paste0("log2FoldChange in ", samplenames[2]))+
      theme_minimal()+
      theme(legend.title = element_blank())
}

compare_significant_genes(res1,res2, samplenames = c("LNCaP", "PC3"), title = "Hypoxia-induced gene expression differences in LNCaP vs PC3 cells")

在这里插入图片描述

我们观察到，确实有一小部分基因是在 LNCaP 和 PC3 比较中均差异表达（绿色标签）。令人满意的是，其中大多数共享基因沿 X = Y 对角线落下。例如，基因，例如 CA9、IGFBP5、STC1 和 PPFIA4 均受到缺氧的严重诱导线。这些基因更有可能是缺氧的真实指标响应。相比之下，仅在一个细胞中具有显著意义的基因线大多位于垂直轴或水平轴上。这些基因是可能与每个细胞背景的特异性反应有关，并且不太可能是缺氧的通用指标。

基因集富集分析
基因集富集分析（GSEA）可用于检测通路是否与另一组相比，由多个基因组成的组中富集。例如，有一个由 200 个基因组成的“基因集” 缺氧途径。我们想确定基因是否与缺氧通路相关的缺氧途径在缺氧条件，作为验证实验是否已执行的一种方式正确。我们还想确定缺氧的其他途径影响细胞内部。

为了进行 GSEA，我们必须根据方式对每个可用的基因进行排名 A组与B组相比，它要高得多或低得多。这称为排名列表，我们已经在上一节中生成了一个使用函数 .我们可以插入这个排名列表从博德研究所（Broad Institute）获得GSEA申请 https://www.gsea-msigdb.org/gsea/index.jsp，或进入 Bioconductor 的 R 包https://bioconductor.org/packages/release/bioc/html/fgsea.html。如果选择使用 GSEA 应用程序而不是 R 包，则排名列表必须采用非常特定的格式。基因必须是所有大写，and 文件应包含一列用于基因名称和具有相应 logfoldChange 的另一列。在这里，我们将使用 R 包执行 GSEA。

可以使用 Bioconductor 安装 R 包命令：

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("fgsea")

我们还需要从 GSEA MSigDB 网页下载路径文件https://www.gsea-msigdb.org/gsea/msigdb/collections.jsp并使用加载它们。我将与以下人员合作 HALLMARK 通路集。
在这里插入图片描述
点击gene symbols 下载，这里我下载的是“h.all.v2023.2.Hs.symbols.gmt”

library(fgsea)
# read in file containing lists of genes for each pathway
hallmark_pathway <- gmtPathways("h.all.v2023.2.Hs.symbols.gmt.txt")
head(names(hallmark_pathway))

[1] "HALLMARK_ADIPOGENESIS"        "HALLMARK_ALLOGRAFT_REJECTION"
[3] "HALLMARK_ANDROGEN_RESPONSE"   "HALLMARK_ANGIOGENESIS"       
[5] "HALLMARK_APICAL_JUNCTION"     "HALLMARK_APICAL_SURFACE"

通路文件作为基因集列表加载。我们可以访问通过使用选择通路名称来获得每个通路的基因：

head(hallmark_pathway$HALLMARK_HYPOXIA, 20)

 [1] "ACKR3"    "ADM"      "ADORA2B"  "AK4"      "AKAP12"  
 [6] "ALDOA"    "ALDOB"    "ALDOC"    "AMPD3"    "ANGPTL4" 
[11] "ANKZF1"   "ANXA2"    "ATF3"     "ATP7A"    "B3GALT6" 
[16] "B4GALNT2" "BCAN"     "BCL2"     "BGN"      "BHLHE40"

加载路径后，我们必须将排名列表格式化为 fgsea 函数可以使用的格式。其中值以基因名称命名。我们还需要处理重复基因名和缺失值。这自定义函数将执行这些操作。

# load the ranked list
#lncap_ranked_list <- read.table("LNCAP_Hypoxia_vs_LNCAP_Normoxia_rank.rnk", header = T, stringsAsFactors = F)
#这里代码会报错“错误于scan(file = file, what = what, sep = sep, quote = quote, dec = dec, : line 184 did not have 2 elements”，以下是我更改的
lncap_ranked_list <- read.delim("LNCAP_Hypoxia_vs_LNCAP_Normoxia_rank.rnk", header = TRUE, stringsAsFactors = FALSE, sep = "\t")
head(lncap_ranked_list)

  Gene.name log2FoldChange
1   CYP4F22       7.550623
2     TGFB1       7.422230
3     GGTA1       7.390595
4     BEND5       7.066944
5      CCN4       6.475747
6     KRT72       6.442253

# formats the ranked list for the fgsea() function
prepare_ranked_list <- function(ranked_list) { 
  # 如果存在重复的基因名，则对其值取平均
  if( sum(duplicated(ranked_list$Gene.name)) > 0) {
    ranked_list <- aggregate(.~Gene.name, FUN = mean, data = ranked_list)
    ranked_list <- ranked_list[order(ranked_list$log2FoldChange, decreasing = T),]#按 log2FoldChange 进行降序排序。
  }
  # 删除包含NA值的行
  ranked_list <- na.omit(ranked_list)
  # 将数据框转换为命名向量
  ranked_list <- tibble::deframe(ranked_list)
  ranked_list
}

lncap_ranked_list <- prepare_ranked_list(lncap_ranked_list)
head(lncap_ranked_list)

 CYP4F22    TGFB1    GGTA1    BEND5     CCN4    KRT72 
7.550623 7.422230 7.390595 7.066944 6.475747 6.442253

现在我们有了命名向量，我们可以将其与标志路径对象一起插入函数。这将生成一个包含与每个相关的浓缩分数的结果表路。

# generate GSEA results table using fgsea() by inputting the pathway list and ranked list
fgsea_results <- fgsea(pathways = hallmark_pathway,
                  stats = lncap_ranked_list,
                  minSize = 15,
                  maxSize = 500,
                  nperm= 1000)
# 排序并选择感兴趣的列
fgsea_results %>% arrange (desc(NES)) %>% select (pathway, padj, NES) %>% head()

                                      pathway       padj      NES
                                       <char>      <num>    <num>
1:                           HALLMARK_HYPOXIA 0.01117818 2.177031
2:                      HALLMARK_ANGIOGENESIS 0.01117818 1.925468
3:                 HALLMARK_ANDROGEN_RESPONSE 0.01117818 1.844683
4:                HALLMARK_TGF_BETA_SIGNALING 0.01117818 1.840881
5:                        HALLMARK_GLYCOLYSIS 0.01117818 1.831385
6: HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION 0.01117818 1.753439

归一化浓缩分数（NES）告诉我们浓缩了多少与常氧样品相比，该途径在缺氧样品中。正如预期的那样，我们观察到最丰富的途径对缺氧是HALLMARK_HYPOXIA途径。

我们可以使用“瀑布”来可视化每个路径的统计数据图，这是的归一化富集分数的侧向条形图每个通路，按显着性进行颜色编码。这个情节很棒用于快速识别显着富集的途径。

waterfall_plot <- function (fsgea_results, graph_title) {
  fgsea_results %>% 
    mutate(short_name = str_split_fixed(pathway, "_",2)[,2])%>% # 移除路径标题中的 'HALLMARK_'
    ggplot( aes(reorder(short_name,NES), NES)) +
      geom_bar(stat= "identity", aes(fill = padj<0.05))+
      #scale_fill_manual(values = c("blue", "red"), name = "padj < 0.05") +
      coord_flip()+
      labs(x = "Hallmark Pathway", y = "Normalized Enrichment Score", title = graph_title)+
      theme(axis.text.y = element_text(size = 7), 
            plot.title = element_text(hjust = 1))
}

waterfall_plot(fgsea_results, "Hallmark pathways altered by hypoxia in LNCaP cells")

在这里插入图片描述
作者原图：

除了缺氧途径外，androgen response和 MTORC1 信号通路似乎在缺氧中也显着富集条件。一些有趣的代谢变化也正在发生。与Warburg效应类似，我们观察到glycolysis 在缺氧条件下大大富集，并且oxidative phosphorylation是负富集的。由于肿瘤变成在癌症进展期间进行性缺氧，该数据可能表明缺氧条件是代谢转换的驱动因素氧化磷酸化为糖酵解作为能量来源。

有趣的是，interferon response pathways都是负面的在缺氧条件下富集，表明缺氧细胞可能对干扰素的反应能力较差。鉴于细胞是培养的在体外（可能不存在干扰素的环境），我发现这个结果出乎意料。

我们还可以通过绘制来缩小感兴趣的特定路径的范围 “富集曲线”使用函数。在在这些图中，X 轴上的黑色刻度表示通路，绿色曲线是衡量基因富集程度的指标对于缺氧组（左侧）或常氧组（右侧）侧）。在缺氧或缺氧中高度富集的通路的示例曲线常氧氧如下图所示。

# 函数封装
plot_enrichment <- function (geneset, pathway, ranked_list) {
  plotEnrichment(geneset[[pathway]], ranked_list)+labs (title = pathway)
}
# 示例的正向富集路径（在缺氧中上调）
plot_enrichment(hallmark_pathway, "HALLMARK_GLYCOLYSIS" , lncap_ranked_list)

在这里插入图片描述

# 示例的反向富集路径（在缺氧中下调）
plot_enrichment(hallmark_pathway, "HALLMARK_OXIDATIVE_PHOSPHORYLATION" , lncap_ranked_list)

在这里插入图片描述

这些类型的富集曲线通常显示在科学中强调上调/下调特定法规的出版物感兴趣的途径。

结论

在这里，我向您展示了一个完整的批量RNA测序分析流程，从下载和处理原始测序文件开始为差异表达基因创建可视化的方法和途径。通过我们的分析，我们观察到：

LNCaP 细胞和 PC3 细胞的全局基因差异很大表达式配置文件。
两种细胞系在缺氧条件下都会上调一组共享的基因条件。
缺氧会上调与糖酵解相关的基因，并且下调与氧化磷酸化相关的基因，提示正在发生代谢转换。

我们已经为每个细胞系生成了差异表达基因列表，可以探索与缺氧相关的新基因响应。此外，我们还学会了如何生产几种不同的可用于批量显示和浏览的可视化效果类型 RNA测序数据。

sepasure

关注

11
点赞
踩
16

收藏

觉得还不错? 一键收藏
0
评论
bs自用教程

参考以下示例：我们的目标是找到差异表达的基因，以响应 LNCaP 和 PC3 细胞系的缺氧。因此，我们将选择两种细胞系的对照样品（Empty_Vector LNCaP 和 siCtrl 的对照样品对于PC3），在常氧或缺氧条件下。具体下表概述了我们需要下载的示例：NCBI以序列读取存档的形式存储测序数据（SRA）文件。每个样本都与一组 SRA 相关联入藏号，如上所示。首先，我们需要下载SRA 针对每个示例运行。然后，我们将使用 SRA 文件生成 FASTQ 文件。
复制链接

扫一扫