Chip-seq上游分析流程学习(一)

这次用到的数据集是GSE274995，里面包含了3个样本的头颈部鳞癌细胞系(Cal27细胞)数据。

数据上传者采用的上游分析软件：

使用了Trimmomatic进行了数据质控；Bowtie2进行了比对分析；SAMtools处理比对后的BAM文件；MACS2用于Peaks Calling；bamCoverage转换为二进制的BigWig文件; 使用了hg19作为参考基因组，现在基本上首先使用GRCh38。

再来看一下曾老师推文和B站视频中的分析流程，其中数据数据过滤环节使用了fastp或者trim-galore软件(替代Trimmomatic)，其他的基本差不多。

总体的步骤可以参考既往mRNA/miRNA的步骤：环境部署——数据下载——查看数据(非质控)——数据质控清洗——数据比对——数据定量，基本逃不出这个框架，只是针对不同类型的数据稍有不同。

本次分析步骤包括：环境部署——数据下载——查看数据(非过滤)

分析步骤

1.环境部署/软件安装：

尝试使用ARM架构(M1/M2芯片)去安装fastqc trim-galore hisat2 subread multiqc samtools salmon fastp，但这些软件中有几个是不兼容的。所以需要改回原来的x86_64架构(Intel芯片)，如果非mac/M1/M2的不需要用这种方式。

# macbook
CONDA_SUBDIR=osx-64 conda create -n chipseq_x86_64 python=3.9
conda activate chipseq_x86_64

# 如果是在服务器中
conda create -n chipseq python=3.9
conda activate chipseq

如果安装软件很慢的话可以先安装mamba再安装软件，比如：只需要把conda修改成mamba即可

conda install -c conda-forge mamba
mamba install -y sra-tools trim-galore samtools deeptools homer meme macs2 bowtie bowtie2 

2.安装软件

conda install -y sra-tools trim-galore samtools deeptools
conda install -y homer meme macs2 bowtie bowtie2

3.下载数据

# 建议在自己的服务器或者本地先建好相应的文件夹
# 然后需要先cd到目标文件夹下
# 不建议直接复制，因为可能会有不知名的空格hhh
nohup cat SRR_Acc_List.txt | while read id; do prefetch $id; done > download.log 2>&1 &

1. nohup的作用是使后续命令在后台运行，即使关闭中断或断开会话任务也会继续进行。
2. cat可以读取文件的内容并将其逐行输出
3. “;”和“|”区别是前者可以分隔独立命令，每条命令按顺序执行，前后无依赖关系。后者是链接命令，前一个命令的输出是后一个命令的输入，命令之间存在数据流依赖。
4. 这里为什么可以用id去作为读取对象，是因为read的命令设计会从输入流中读取每一行的数据，并存储到指定的变量中，id是自定义的变量名，假设换成line或者其他都可以。
5. ">"的作用是输入重定向，也就是说如果在download.log前面使用了比如“|”，这个log文件会被当做命令。
6. “2>&1”的含义是2(表示标准错误，stderr，文件描述2)，>(表示重定向)，&1(表示重定向的目标是标准输出，stdout, 文件描述1)。
7. 最后的“&”是将整个命令放在后台运行。

4.sra文件转化为fastq

先看看样本的单双链信息，点击具体样本

从这里可以看到样本采用的测序方式是什么, Layout: SINGLE,说明是单端测序数据

此外在这里稍提一个概念，有时候我们在让公司帮忙做测序的时候，公司人员交流会提到“测多少个G”这个概念。这里的测多少个G是表示Gigabases(表示Giga碱基，10亿碱基)。

我们可以看到下面这个页面上最上面的这句话：1 ILLUMINA (Illumina NovaSeq 6000) run: 42.3M spots, 3.2G bases, 1.2Gb downloads。

1. ILLUMINA (Illumina NovaSeq 6000): ILLUMINA：指使用的测序技术平台，Illumina 是一种高通量测序技术。Illumina NovaSeq 6000：是 Illumina 平台的一种具体型号。
2. run:表示这是一次测序实验（run）的统计数据。
3. 42.3M spots:42.3M：表示测序的 reads（读取数），单位是百万（M，即 10^6）。spots：指的是测序数据中的原始数据点，通常等同于测序读取数（reads），或者是指成对读取的原始数据。
4. 3.2G bases:3.2G：表示测序生成的碱基数，总量是 3.2 亿个碱基（G = 10^9）。bases：指测序生成的碱基序列。这里G bases的简写也是 Gb，要与下面的Gigabytes作区分，这里表示碱基数。
5. 1.2Gb downloads:1.2Gb：表示这次测序运行生成的文件大小约为 1.2 吉字节（Gb = gigabytes）。downloads：指该测序运行的数据文件总大小为 1.2GB，可以通过 SRA 下载，这里的Gb表示存储量大小。
6. 一般来说，测序数目(42.310^6) * 测序长度(?这里是未知) = 测序量(3.210^9 bases)，我们反推一下测序长度，当然这个数据是单端测序数据，如果是双端我们还需要除以2。

我们用一种不太智能的方法验证一下是不是75bp，我们进入如下界面，点开reads。

数一数确实是75个碱基。

接下来正式开始进行fastq转换

# 由于笔者这里存储的有点乱，所以会比较复杂
# 使用者在自行使用的时候，如果觉得设置路径很麻烦就把文件全放一个文件夹下面，先学会再说
# 定义变量，指定工作目录路径
# path路径需要通过pwd确定

# 可以选择写一个shell脚本
# 把下面代码输入进去
#!/bin/bash

# 也可以直接复制这些代码进行运行
path="/home/lm/Z_Projects/chipseq"
list_file="${path}/SRR_Acc_List.txt"  # 列表文件路径
sra_dir="${path}/sra/sra"             # 输入的 .sra 文件所在路径
fastq_dir="${path}/fastq"             # 输出的 .fastq 文件存储路径

# 创建 fastq 输出目录
mkdir -p "${fastq_dir}"

# 读取列表文件并处理
cat "${list_file}" | while read id; do
  # 检查 .sra 文件是否存在
  sra_file="${sra_dir}/${id}.sra"
  if [ -f "${sra_file}" ]; then
    echo "Processing ${id}..."
    # 使用 fastq-dump 处理 .sra 文件，输出到 fastq 文件夹
    fastq-dump --gzip -O "${fastq_dir}" "${sra_file}"
    echo "Completed ${id}."
  else
    echo "File ${sra_file} not found, skipping..."
  fi
done

1. 读取列表文件并循环处理：cat "读取{list_file} 变量指定的列表文件内容，通常是一个包含若干 SRR ID（例如 SRR123456）的文本文件，每一行是一个 ID。|: 管道符，将 cat 的输出作为后续 while 循环的输入。while read id: 按行读取列表文件中的内容，并将每一行存储到变量 id 中，逐一处理。do: 标志循环的开始，后续是循环体的内容。
2. 定义 .sra 文件的路径：定义文件所在的目录路径。{id}: 当前循环中读取的 ID 值，例如 SRR123456。"{id}.sra": 拼接出完整的 .sra 文件路径，例如 /home/user/chipseq/sra/sra/SRR123456.sra。
3. 检查 .sra 文件是否存在：[ -f "检查文件是否存在且是一个普通文件。条件判断语句，如果文件存在，则执行之后的代码块；否则跳过。{sra_file}": 表示要检查的文件路径。
4. 打印处理消息：echo: 打印消息到终端，告知当前正在处理的 ID。
5. 使用 fastq-dump 工具处理文件：fastq-dump: 是一个常用的 SRA 数据转换工具，用于将 .sra 格式转换为 .fastq 或 .fastq.gz 格式。--gzip: 转换后的 .fastq 文件会直接压缩为 .fastq.gz 格式，节省存储空间。-O "指定输出目录，将转换后的文件存储到{fastq_dir} 目录中。"${sra_file}": 要处理的 .sra 文件路径。
6. 打印完成消息：echo: 打印消息到终端，告知用户当前 ID 的处理已经完成。文件不存在时的提示：else: 如果条件 if 不成立（即文件不存在），执行该分支的代码。echo: 打印警告消息到终端，告知用户某个 ID 的 .sra 文件未找到，并跳过该文件的处理。
7. 循环结束标志：fi: 结束 if 条件判断块。done: 表示 while 循环的结束。

4.查看数据(非过滤)

# cd到fastq文件夹中

# 不整合，输出每一个样本的质控报告
# fastqc -t [线程数] [存储路径] [文件来源路径]
fastqc -t 6 -o ./ ./SRR*.fastq.gz

不整合，输出每一个样本的质控报告

fastqc运行 FastQC 软件; -t 6 运行6个线程，加快速度；./ 存在当前文件夹中； ./SRR*.fastq.gz 读取fastq.gz文件

可以打开html文件看一看，具体的质控报告细节可以看一下既往的转录组上游分析流程~ 也可以看看其他UP主的推文。

# 整合，输出一个总的质控报告
# 请注意，这里是先要进行单一质控，然后根据单一质控得到的zip文件结果进行整合！
multiqc  ./*.zip  -o ./

整合，输出一个总的质控报告。这里是先要进行单一质控，然后根据单一质控得到的zip文件结果进行整合！

既往推文

1. 转录组上游分析流程(一）：https://mp.weixin.qq.com/s/bwUFJ-kBdUTp9WyQRQPnyg
2. 转录组上游分析流程(二）：https://mp.weixin.qq.com/s/Lq9lES4M6ZzzleXfg0iWCw
3. 转录组上游分析流程(三）：https://mp.weixin.qq.com/s/EIoPqIKC4IuOFChCmHlsVQ
4. 转录组上游分析流程(四）：https://mp.weixin.qq.com/s/jYGltNVfdRu-k28Yhxlllg

参考资料

1. 生信技能树：https://mp.weixin.qq.com/s/KOFd60USQdKY2C4Sc4OIvA
2. 生信学习日记：https://mp.weixin.qq.com/s/z-qbVZqo-cKT3Tfv4LT6OQ
3. Tobin笔记：https://mp.weixin.qq.com/s/Uf8GgK4kZAWbRKKoh6IiyA

致谢：感谢曾老师以及生信技能树团队全体成员。

注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟

- END -