咱们《生信技能树》的B站有一个ChIP-seq数据分析实战视频课程,缺乏配套笔记。恰好前些天的求贤令:曾经我给你带来了十万用户,但现在祝你倒闭,走大运结识了几位优秀小伙伴!
其中中国医科大的“小高”同学给大家带来的就是ChIP-seq数据分析实战视频课程的配套笔记,希望可以帮助大家更好的吸收消化课程内容!
首先视频免费共享在B站:【生信技能树】Chip-seq测序数据分析ChIP-SEQ实战演练的素材:链接:https://share.weiyun.com/53CwQ8B 密码:ju3rrh, 包括一些公司PPT,综述以及文献 ChIP-SEQ 实战演练的思维导图:文档链接:https://mubu.com/doc/11taEb9ZYg 密码:wk29
接下来是根据生信技能树课程、思维导图、PPT和综述文献整理的笔记。前面的笔记是:ChIP-seq数据分析课程学习笔记之fastq测序数据的准备
前面我们已经下载好了全部的fastq数据,接下来就
先走一波QC
cd /data/gy/project/epi/qc
## 相对目录需要理解
ls ../raw/*gz|xargs fastqc -t 10 -o ./
官网上有个相关的视频教程:Using FastQC to check the quality of high throughput sequence,链接:https://www.youtube.com/watch?v=bz93ReOv87Y。(我上传了百度云链接: https://pan.baidu.com/s/1El3ZCli9uVC7u-58tj3HBw 密码: 2wj2)帮助文档在:https://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/。
左侧可以看到大致的结果:绿色PASS、黄色Warning、红色Failure
包含基本信息:测序平台的版本和相应的编码版本号、reads的数量、长度、GC含量等信息。
所有碱基的测序质量的统计结果。
此图横轴是测序序列1-36个碱基。纵轴是质量得分,20表示0.01的错误率,30表示0.001。
箱线图中红线表示中位数,蓝线是平均数的连线。
Warning 报警:任一碱基质量<10,或者是任一中位数<25
Failure 报错:任一碱基质量<5,或者是任一中位数<20
各个tile的序列质量情况。此图横轴是测序序列1-36个碱基。纵轴是tile的编号。
蓝色表示测序质量很好,颜色越红越不好,后续分析可以把这个tile结果去除。
各个序列质量的分布。横图表示Q值(测序质量quality),纵轴是read数目。绝大多数序列质量在30以上可以说质量很好。
每条序列中各个位置的平均碱基比例,如出现AT或GC分离的情况说明这个数据有问题,需要处理。
序列平均GC含量分布图。横轴是GC比例0 - 100%;纵轴是每条序列GC含量对应的数量。
蓝色:理论值。红色:真实值。两条越接近约好。
序列中各个位点的N含量(测序机器识别不出到底是何种碱基时,会产生N),越小越好。
序列测序长度统计。每次测序仪测出来的长度在理论上应该是完全相等的,但是实际总会有一些偏差。此图中,36bp是是主要的,但是还是有少量的35和37bp。当测序的长度严重不同时,表明此次测序不成功。
统计序列的重复度。低重复水平可能表明目标序列的覆盖率非常高,但高重复水平可能表明存在某种富集偏倚,例如 PCR over amplification。
某个序列大量出现是over-represented。
这个图表示序列中两端adapter的情况。一般测序仪下机之前会去除接头序列。
友情提示MultiQC安装可能会遇到坑:
依赖python2.7+, 3.4+ 或者 3.5+
#pip安装
pip install git+https://github.com/ewels/MultiQC.git
# conda安装
conda install -y bioconda multiqc
直接指定MultiQC要分析的文件路径即可,若数据在当前目录下输入multiqc ./
即可。
以上为简要介绍,更多精彩内容可以参加**生信入门课**
这个时候选择trim_galore软件进行过滤,过滤条件:测序得到的原始序列含有接头序列或低质量序列,为了保证信息分析的准确性, 需要对原始数据进行质量控制,得到高质量序列(即Clean Reads),原始序列质量控制的标准为:
①去除含接头的reads;
②过滤去除低质量值数据,确保数据质量;
③去除含有N(无法确定碱基信息)的比例大于5%的reads;
单端测序数据的代码如下;
cd /data/gy/project/epi/clean
#改成自己的
analysis_dir=/data/gy/project/epi
bin_trim_galore="trim_galore"
ls ../raw/*gz | while read fq1;
do
nohup $bin_trim_galore -q 25 --phred33 --length 25 -e 0.1 --stringency 4 -o $analysis_dir/clean $fq1 &
done
#再QC一下
ls ../clean/*gz|xargs fastqc -t 10 -o ./
multiqc ./
然后直接用bowtie2进行比对和统计比对率, 需要提前下载参考基因组然后使用命令构建索引,或者直接就下载索引文件:
下载小鼠参考基因组的索引和注释文件, 这里用常用的mm10
# 索引大小为3.2GB, 不建议自己下载基因组构建,可以直接下载索引文件,代码如下:
mkdir referece && cd reference
wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip
unzip mm10.zip
单端测序数据的比对代码如下:
cd /data/gy/project/epi/align
## 相对目录需要理解
bin_bowtie2='/root/miniconda3/envs/chipseq/bin/bowtie2'
bin_bowtie2=bowtie2
bowtie2_index=/data/gy/public/reference/index/bowtie/mm10
## 一定要搞清楚自己的bowtie2软件安装在哪里,以及自己的索引文件在什么地方!!!
ls ../clean/*gz |while read id;
do
file=$(basename $id )
sample=${file%%.*}
echo $file $sample
$bin_bowtie2 -p 5 -x $bowtie2_index -U $id | samtools sort -O bam -@ 5 -o - > ${sample}.bam
done
对bam文件进行QC
cd /data/gy/project/epi/align
ls *.bam |xargs -i samtools index {}
ls *.bam | while read id ;do (nohup samtools flagstat $id > $(basename $id ".bam").stat & );done
比对成功率都挺好的:
Control_1_trimmed.stat:7438540 + 0 mapped (88.03% : N/A)
Control_2_trimmed.stat:7221781 + 0 mapped (86.40% : N/A)
H2Aub1_1_trimmed.stat:8969578 + 0 mapped (97.40% : N/A)
H2Aub1_2_trimmed.stat:13229916 + 0 mapped (97.53% : N/A)
H3K36me3_1_trimmed.stat:11737310 + 0 mapped (98.89% : N/A)
Ring1B_1_trimmed.stat:4634240 + 0 mapped (93.59% : N/A)
Ring1B_2_trimmed.stat:4646919 + 0 mapped (93.85% : N/A)
RNAPII_S2P_1_trimmed.stat:25018794 + 0 mapped (97.26% : N/A)
RNAPII_S2P_2_trimmed.stat:6112834 + 0 mapped (95.00% : N/A)
RNAPII_S2P_3_trimmed.stat:8675514 + 0 mapped (96.99% : N/A)
RNAPII_S5P_2_trimmed.stat:12182274 + 0 mapped (98.17% : N/A)
RNAPII_S5PRepeat_1_trimmed.stat:4163763 + 0 mapped (82.81% : N/A)
RNAPII_S7P_1_trimmed.stat:6386269 + 0 mapped (80.90% : N/A)
RNAPII_S7P_2_trimmed.stat:5971178 + 0 mapped (82.66% : N/A)
题外话:真的是 大神一句话,菜鸟跑半年。这两节加起来20分钟的课程,感觉自己踩的坑大概比这篇笔记的字数还多。经验教训:多读文献,多做实战,shell脚本还要好好看。
联系客服