免费视频课程-表观调控组学联合分析实战

我这七年写了几万篇教程，制作了几百个小时的教学实战演练视频课程，都是免费分享在各大网站（B站，知乎，简书，博客，GitHub，微云），必然会出现部分教程过时，一些资料缺失（主要是链接失效）。而且很多平台都是生信技能树的各个志愿者帮忙打理，我不可能要求志愿者们在辛辛苦苦帮我整理和发布资料的同时还提供答疑。

今天，我们带来的是表观调控组学联合分析实战-B站免费课程交流群，但是并不想为它设置交流群了。

我们所有的群都是由生信技能树的官方拉群小助手操作，大家在应该是多次看到了：

也成功的组建了分门别类的学习交流群，如下：

其实所谓的表观调控组学联合分析实战就是上面已经建立了的《RNA-seq数据分析》搭配一个《ATAC-seq数据分析》或者《ChIP-seq数据分析》。搞那么多群，也确实维护不过来。

课程获取方式

视频首先可以在B站看到：https://www.bilibili.com/video/av77098514 （失效日期，随时），请尽快学习
视频课程配套课件（幕布思维导图）：https://mubu.com/doc/3Bd4aieYug （理论上不会失效）
配套代码目前也是微云下载：https://share.weiyun.com/5qmxu7Z（失效日期，随时），请尽快保存

课程目录

01 果蝇参考基因组和注释文件准备
02 文献测序原始数据下载
03 ChIP-seq 数据质量控制与过滤
04 ChIP-seq 数据从比对到 Peaks calling
05 RNA-seq 数据从比对到定量
06 用幕布展示此课程目录
07 验证目标基因的指定外显子 knock out 是否有效
08 多个样本基因信号矩阵来相关性计算与可视化
09 对RNA-seq数据找差异基因以及可视化
10 根据差异倍数来计算样本间相关性
11 使用 ChIPseeker 对单组 Peak 进行注释与可视化
12 对多组 Peaks 进行批量注释与可视化
13 通过 bw 文件来计算 ChIP-seq 样本间相关性
14 探索ChIPQC 包用法并且查看 ChIP-seq 数据质量
15 最新版参考基因组注释peaks之果蝇 dm6
16 样本间 Peaks 交集韦恩图
17 基因组版本注释转换
18 使用 R 和 linux 对多个 bed 进行骚操作处理
19 对单个 bed 文件和 bw 文件可视化
20 多个 bw 文件进行可视化
21 获得差异表达基因 bed 文件
22 对多个基因集进行富集分析

需要生信入门基础知识

当然需要读者具备计算机基础知识，我把它粗略的分成基于R语言的统计可视化，以及基于Linux的NGS数据处理：

其中，对甲基化芯片数据处理实战来说，R更重要一点，我把R的知识点路线图搞定，如下：

了解常量和变量概念
加减乘除等运算（计算器）
多种数据类型（数值，字符，逻辑，因子）
多种数据结构（向量，矩阵，数组，数据框，列表）
文件读取和写出
简单统计可视化
无限量函数学习

有一些需要更新

生信技能树绝大部分课程，都是三五年前录制的，所以很多软件更新了，那个时候也不怎么流行aspera高速下载！

主要是上游流程的更新。

环境搭建

如果是全新服务器或者全新用户，首先需要安装conda（最适合初学者的软件管理解决方案）：

#一路yes下去
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-4.6.14-Linux-x86_64.sh
source ~/.bashrc

然后使用conda安装一些软件或者软件环境，比如下载测序数据文件的aspera软件环境：

conda create -n download -y
conda activate download
conda install -y -c hcc aspera-cli
which ascp
## 一定要搞清楚你的软件被conda安装在哪
ls -lh ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh

比如ATAC-SEQ数据分析流程的相关软件：

## 安装好conda后需要设置镜像。
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes

conda create -n atac -y python=2 bwa
conda info --envs
# 可以用search先进行检索
conda search trim_galore
## 保证所有的软件都是安装在 atac 这个环境下面
conda activate atac
conda install -y trim-galore bedtools deeptools homer meme macs2 bowtie bowtie2 sambamba
conda search samtools
conda install -y samtools=1.11

然后构建工作目录架构：

# 注意组织好自己的项目
mkdir -p ~/project/atac/
cd ~/project/atac/
mkdir {sra,raw,clean,align,peaks,motif,qc}
cd raw

取决于个人习惯。

实战数据准备

参考：使用ebi数据库直接下载fastq测序数据 , 需要自行配置好，然后去EBI里面搜索到的 fq.txt 路径文件：

项目地址是

#一次性下载所有的 fastq.gz样本
dsa=$HOME/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh
ls -lh $dsa
# conda activate download
# 自己搭建好 download 这个 conda 的小环境哦。
x=_1
y=_2
for id in {73..80}
do
 ascp -QT -l 300m -P33001 -i $dsa \
 era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR123/0$id/SRR123031$id/SRR123031$id$x.fastq.gz . 
 ascp -QT -l 300m -P33001 -i $dsa \
 era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR123/0$id/SRR123031$id/SRR123031$id$y.fastq.gz .
done

把上面的代码存为代码文件： download.sh ，然后使用下面的命令放在后台下载即可：

conda activate download
nohup bash download.sh &

得到的文件如下：

1.2G 10月 9 22:07 SRR12303173_1.fastq.gz
1.2G 10月 9 22:12 SRR12303173_2.fastq.gz
1.1G 10月 9 22:17 SRR12303174_1.fastq.gz
1.2G 10月 9 22:24 SRR12303175_1.fastq.gz
1.3G 10月 9 22:30 SRR12303175_2.fastq.gz
1.1G 10月 10 09:54 SRR12303176_1.fastq.gz
1.1G 10月 10 09:56 SRR12303176_2.fastq.gz
1.4G 10月 10 10:02 SRR12303177_1.fastq.gz
1.4G 10月 10 10:05 SRR12303177_2.fastq.gz
1.2G 10月 10 10:09 SRR12303178_1.fastq.gz
1.2G 10月 10 10:13 SRR12303178_2.fastq.gz
1.2G 10月 10 10:18 SRR12303179_1.fastq.gz
1.2G 10月 10 10:21 SRR12303179_2.fastq.gz
1.3G 10月 10 10:26 SRR12303180_1.fastq.gz
1.3G 10月 10 10:35 SRR12303180_2.fastq.gz

可以看到，aspera下载的时候，中间11个小时任务终止了，是我自己重新跑了aspera下载，续起来了的。而且如果你仔细看会发现 SRR12303174这个样品只有R1的fq文件缺失了R2，也是需要重新单独下载。

conda activate download
dsa=$HOME/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh
id=74;y=_2
ascp -QT -l 300m -P33001 -i $dsa \
era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR123/0$id/SRR123031$id/SRR123031$id$y.fastq.gz .
## SRR12303174_2.fastq.gz 100% 1094MB 87.3Mb/s 01:27
## Completed: 1120670K bytes transferred in 87 seconds
## (104410K bits/sec), in 1 file.
## 全部文件下载完毕后，使用下面的命令检查一下fq.gz文件是否完整
gzip -t *gz

只有项目的fq数据全部准备而且确认无误后才能进行下一步！

测序数据的质量控制

这里选择trim_galore软件，自动批量运行：

mkdir -p ~/project/atac/
cd ~/project/atac/
conda activate atac
trim_galore --help
for id in {73..80}
do
 nohup trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired \
 -o clean raw/SRR123031$id*.fastq.gz &
done

得到的文件如下：

 828M 10月 10 17:00 clean/SRR12303173_1_val_1.fq.gz
 836M 10月 10 17:00 clean/SRR12303173_2_val_2.fq.gz
 797M 10月 10 18:09 clean/SRR12303174_1_val_1.fq.gz
 808M 10月 10 18:09 clean/SRR12303174_2_val_2.fq.gz
 900M 10月 10 17:03 clean/SRR12303175_1_val_1.fq.gz
 917M 10月 10 17:03 clean/SRR12303175_2_val_2.fq.gz
 787M 10月 10 16:58 clean/SRR12303176_1_val_1.fq.gz
 794M 10月 10 16:58 clean/SRR12303176_2_val_2.fq.gz
 992M 10月 10 17:13 clean/SRR12303177_1_val_1.fq.gz
1006M 10月 10 17:13 clean/SRR12303177_2_val_2.fq.gz
 845M 10月 10 17:01 clean/SRR12303178_1_val_1.fq.gz
 855M 10月 10 17:01 clean/SRR12303178_2_val_2.fq.gz
 840M 10月 10 17:01 clean/SRR12303179_1_val_1.fq.gz
 847M 10月 10 17:01 clean/SRR12303179_2_val_2.fq.gz
 945M 10月 10 17:06 clean/SRR12303180_1_val_1.fq.gz
 963M 10月 10 17:06 clean/SRR12303180_2_val_2.fq.gz

这个过滤还是有点狠的，之前1.3G现在都小于1G了。实际上可以走fastqc+multiqc的质控看过滤前后的具体情况。

数据比对到参考基因组

1、mm10小鼠参考基因组的下载

#下载
mkdir -p ~/project/atac/ref
cd ~/project/atac/ref
nohup wget http://hgdownload.soe.ucsc.edu/goldenPath/mm10/bigZips/mm10.fa.gz &
#解压
gunzip mm10.fa.gz

2、bowtie2-build构建参考基因组索引文件

这一步会生成6个索引文件,这一步耗时比较常。也可以自行下载对应的参考基因组索引。

conda activate atac
nohup bowtie2-build mm10.fa mm10 1>log 2>&1 &

得到的文件如下：

848M 10月 10 17:20 mm10.1.bt2
633M 10月 10 17:20 mm10.2.bt2
6.0K 10月 10 16:36 mm10.3.bt2
633M 10月 10 16:36 mm10.4.bt2
2.6G 1月 23 2020 mm10.fa
848M 10月 10 18:05 mm10.rev.1.bt2
633M 10月 10 18:05 mm10.rev.2.bt2

3、bowtie2进行批量比对

首先制作配置文件：

cd ~/project/atac/align

ls $HOME/project/atac/clean/*_1.fq.gz > 1
ls $HOME/project/atac/clean/*_2.fq.gz > 2
ls $HOME/project/atac/clean/*_2.fq.gz |cut -d"/" -f 7|cut -d"_" -f 1 > 0
paste 0 1 2 > config.clean ## 供mapping使用的配置文件

然后创建含有如下内容的脚本（aligh.sh）：

## 相对目录需要理解
bowtie2_index=$HOME/project/atac/ref/mm10
## 一定要搞清楚自己的bowtie2软件安装在哪里，以及自己的索引文件在什么地方！！！
#source activate atac
cat config.clean |while read id;
do echo $id
arr=($id)
fq2=${arr[2]}
fq1=${arr[1]}
sample=${arr[0]}
## 比对过程15分钟一个样本
bowtie2 -p 5 --very-sensitive -X 2000 -x $bowtie2_index -1 $fq1 -2 $fq2 |samtools sort -O bam -@ 5 -o - > ${sample}.raw.bam
samtools index ${sample}.raw.bam
bedtools bamtobed -i ${sample}.raw.bam > ${sample}.raw.bed
samtools flagstat ${sample}.raw.bam > ${sample}.raw.stat
# https://github.com/biod/sambamba/issues/177
sambamba markdup --overflow-list-size 600000 --tmpdir='./' -r ${sample}.raw.bam ${sample}.rmdup.bam
samtools index ${sample}.rmdup.bam

## ref:https://www.biostars.org/p/170294/
## Calculate %mtDNA:
mtReads=$(samtools idxstats ${sample}.rmdup.bam | grep 'chrM' | cut -f 3)
totalReads=$(samtools idxstats ${sample}.rmdup.bam | awk '{SUM += $3} END {print SUM}')
echo '==> mtDNA Content:' $(bc <<< "scale=2;100*$mtReads/$totalReads")'%'

samtools flagstat ${sample}.rmdup.bam > ${sample}.rmdup.stat
samtools view -h -f 2 -q 30 ${sample}.rmdup.bam |grep -v chrM |samtools sort -O bam -@ 5 -o - > ${sample}.last.bam
samtools index ${sample}.last.bam
samtools flagstat ${sample}.last.bam > ${sample}.last.stat
bedtools bamtobed -i ${sample}.last.bam > ${sample}.bed
done

提交脚本的代码是：

conda activate atac
nohup bash aligh.sh 1>log 2>&1 &

全部运行完毕后输出非常多文件。

首先看bam文件，如下：

1.1G 10月 11 15:49 SRR12303173.last.bam
1.8G 10月 10 23:01 SRR12303173.raw.bam
1.3G 10月 11 15:48 SRR12303173.rmdup.bam
823M 10月 11 16:00 SRR12303174.last.bam
1.7G 10月 11 00:24 SRR12303174.raw.bam
976M 10月 11 15:59 SRR12303174.rmdup.bam
1.4G 10月 11 16:11 SRR12303175.last.bam
2.2G 10月 11 02:26 SRR12303175.raw.bam
1.6G 10月 11 16:09 SRR12303175.rmdup.bam
1.2G 10月 11 16:23 SRR12303176.last.bam
1.8G 10月 11 03:52 SRR12303176.raw.bam
1.4G 10月 11 16:21 SRR12303176.rmdup.bam
1.8G 10月 11 16:37 SRR12303177.last.bam
2.5G 10月 11 06:16 SRR12303177.raw.bam
2.1G 10月 11 16:35 SRR12303177.rmdup.bam
1.2G 10月 11 16:50 SRR12303178.last.bam
1.9G 10月 11 07:53 SRR12303178.raw.bam
1.4G 10月 11 16:48 SRR12303178.rmdup.bam
1.2G 10月 11 17:02 SRR12303179.last.bam
1.9G 10月 11 09:35 SRR12303179.raw.bam
1.4G 10月 11 17:00 SRR12303179.rmdup.bam
1.7G 10月 11 17:16 SRR12303180.last.bam
2.4G 10月 11 11:51 SRR12303180.raw.bam
1.9G 10月 11 17:14 SRR12303180.rmdup.bam

每个样品分别会输出3个bam文件，测序数据比对的bam，以及去除PCR重复后的bam，以及去除线粒体reads后的bam文件。

查看log日志，发现这些样本的线粒体含量是：

==> mtDNA Content: 1.81%
==> mtDNA Content: 3.72%
==> mtDNA Content: 1.88%
==> mtDNA Content: 1.98%
==> mtDNA Content: 2.16%
==> mtDNA Content: 3.78%
==> mtDNA Content: 2.11%
==> mtDNA Content: 2.17%

因为我们是首先去除PCR重复然后计算线粒体含量，其实是不准确的。

比对后的bam文件的统计

测序文库复杂度的检验

一个简单的含有awk脚本的shell脚本，代码如下：

ls *.last.bam|while read id;
do
 bedtools bamtobed -bedpe -i $id | \
 awk 'BEGIN{OFS="\t"}{print $1,$2,$4,$6,$9,$10}' | sort | uniq -c | \
 awk 'BEGIN{mt=0;m0=0;m1=0;m2=0} ($1==1){m1=m1+1} ($1==2){m2=m2+1} {m0=m0+1} {mt=mt+$1} END{m1_m2=-1.0; if(m2>0) m1_m2=m1/m2;printf "%d\t%d\t%d\t%d\t%f\t%f\t%f\n",mt,m0,m1,m2,m0/mt,m1/m0,m1_m2}' > ${id%%.*}.nodup.pbc.qc;
done

脚本制作好了后命名为：

conda activate atac
nohup bash stat_qc.sh &

Library complexity measures计算结果如下，…nodup.pbc.qc文件格式为：

TotalReadPairs
DistinctReadPairs
OneReadPair
TwoReadPairs
NRF=Distinct/Total
PBC1=OnePair/Distinct
PBC2=OnePair/TwoPair

针对NRF、PBC1、PBC2这几个指标，ENCODE官网提供了标准.

计算结果显示NRF、PBC1、PBC2的值都非常完美，说明我们进行过滤和PCR去重的bam文件质量上没有问题，可以用于后续的分析。

生信菜鸟团

欢迎去论坛biotrainee.com留言参与讨论，或者关注同名微信公众号biotrainee

免费视频课程-表观调控组学联合分析实战

课程获取方式

课程目录

需要生信入门基础知识

有一些需要更新

环境搭建

实战数据准备

测序数据的质量控制

数据比对到参考基因组

1、mm10小鼠参考基因组的下载

2、bowtie2-build构建参考基因组索引文件

3、bowtie2进行批量比对

比对后的bam文件的统计

测序文库复杂度的检验

2024年11月
一	二	三	四	五	六	日
« 九
				1	2	3
4	5	6	7	8	9	10
11	12	13	14	15	16	17
18	19	20	21	22	23	24
25	26	27	28	29	30