写在前面的话
新学期学新东西，接触表观遗传方面，才发现对基因组，转录组的知识十分欠缺，留下学习笔记也供自己复习。文末有联系方式，最近有不少私信，直接挂出来好了！！！欢迎交流！
以m6A-seq为例，梳理相关分析流程。
表观转录组简介
表观转录组：主要研究RNA所携带的化学修饰对基因表达的影响，简而言之就是针对AGCU碱基上发生的修饰碱基类型。
目前发现超过160种修饰类型标红表示目前有成熟的检测技术
mRNA上的RNA修饰
引自：Li, X., Xiong, X. & Yi, C. Epitranscriptome sequencing technologies: decoding RNA modifications. Nat Methods 14, 23–31 (2017).
例如：m6A：表示在A这个碱基（腺嘌呤）的第六个的第6位氮原子（N）发生了甲基化修饰，主要存在终止密码子附近。
典型的蛋白编码基因结构
引自：https://en.wikipedia.org/wiki/File:Gene_structure_eukaryote_2_annotated.svg
高通量表观转录组学测序
包括抗体富集，自身修饰错配特性，酶活性等原理。以m6A举例：
m6A对基因转录表达的动态调控包括：RNA在转录过程中就会加上m6A修饰；m6A修饰调控可变剪切；m6A修饰调控翻译
研究文章类型：1.表观转录组机制研究（鉴定新的表观转录机器）；2.新的RNA修饰检测方法（新的修饰/新的检测技术/高通量测序技术）；3.不同体系下RNA修饰的功能研究 ；4.不同体系/状态下的修饰图谱绘制（resources） 5.RNA修饰的相关算法（问题导向/数据导向）；6.RNA修饰在临床上的研究 ；7.RNA修饰的数据库构建
m6A-seq测序流程图片段化后，用m6A抗体富集m6A修饰片段，得到两个文库，input也为RNA-seq，IP为m6A-seq。
富集m6A修饰的peak主要位于TSS和Stop codon处。
引自：Dominissini, D., Moshitch-Moshkovitz, S., Schwartz, S. et al. Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature 485, 201–206 (2012).
m6A-seq分析流程
软件安装
conda active m6A-seq
conda install -c bioconda sra-tools # sra-tools用于下载数据
conda install -c bioconda fastqc # fastqc查看测序数据质量
conda install -c bioconda trim_galore # trim_galore去除接头和低质量的 reads
conda install -c bioconda hisat2 # hisat2比对软件
conda install -c bioconda samtools # samtools 处理sam和bam文件
conda install -c bioconda bedtools # bedtools 处理 bed 文件
conda install -c bioconda macs2 # macs2 用于Peak Calling
conda install -c bioconda homer # HOMER注释peak
conda install -c bioconda deeptools # chip-seq综合数据处理
测序原始数据质控
vim 01.trim.sh
for file *.fq.gz
do
trim_galore -q 20 --phred33 --stringency 3 --length 35 -o ./ ${file}
done
sh 01.trim.sh
#trim_galore是调用cutadapt & fastQC的整合质控去接头工具
常用参数
--adapter：输入adapter序列。也可以不输入，Trim Galore会自动寻找可能性最高的平台对应的adapter。
-q reads质量的cutoff
--stringency：设定可以忍受的前后adapter重叠的碱基数，默认为1（非常苛刻）可以放宽
--length：设定输出reads长度阈值，小于设定值会被抛弃。
--paired：对于双端测序结果，一对reads中，如果有一个被剔除，那么另一个会被同样抛弃，而不管是否
达到标准。
--retain_unpaired：对于双端测序结果，一对reads中，如果一个read达到标准，但是对应的另一个要被抛
弃，达到标准的read会被单独保存为一个文件。
--gzip和--dont_gzip：清洗后的数据zip打包或者不打包。
--output_dir：输出目录。需要提前建立目录，否则运行会报错。
--fastqc：处理完成了运行fastQC质控
#nohup sh 02.remove_rRNA.sh > remove_rRNA.sh.log 2>&1 &
2>&1的运用 1：标准输出；2：标准错误。2>&1表明将文件描述2（标准错误输出）的内容重定向到文件描述符1（标准输出）。可以很好的将错误信息保存，帮助我们定位问题。
去除rRNA影响
因IP富集的文库含有较多的rRNA，因rRNA也具有polyA尾，去除rRNA上的m6A修饰影响
vim 02.remove_rRNA.sh
for file in *trimmed.fq.gz
do
hisat2 -p 10 -k 1 -x /human_rRNA_index/rRNA_seq_fix -U ${file} -S ${file}.ht2-rRNA.sam --un-gz ${file}_de-rRNA.fq.gz > ${file}.hisat2.log
done
sh 02.remove_rRNA.sh # 去除rRNA影响
# 构建rRNA序列索引，NCBI下载
基因组回帖比对
vim 03.map_to_genome.sh
for file in *fq.gz_de-rRNA.fq.gz
do
hisat2 -p 10 -k 1 -x /reference_geneme/hg19/hg19_only_chromosome -U ${file} -S ${file}.sam
done
sh 03.map_to_genome.sh # 基因组回帖
# 构建比对基因组索引
hisat2-build $genome_seq.fa $genome_seq
# 参考基因组 & 注释文件下载 Ensembl，UCSC，NCBI
排序转化bam文件
vim 04.samtools_sort.sh
for file in *_de-rRNA.fq.gz.sam
do
samtools sort -@ 5 -O BAM -o ${file}.sort.bam ${file}
done
sh 04.samtools_sort.sh # 排序，转化bam文件
bam文件建立index
vim 05.samtools_index.sh
for file in *de-rRNA.fq.gz.sam.sort.bam
do
samtools index ${file}
done
sh 05.samtools_index.sh # 构建bam文件索引
samtools介绍：SAM文件转换为BAM文件、排序：SAMtools。-@ 线程数 -O 输出文件格式 （BAM / SAM / CRBAM） -o 输出文件 view: BAM-SAM/SAM-BAM 转换和提取部分比对 sort: 比对排序 merge: 聚合多个排序比对 index: 索引排序比对 mpileup: 产生基于位置的结果和 consensus/indel calling
bigwig文件在IGV可视化
vim 06.bam_to_bigwig.sh
for file in *de-rRNA.fq.gz.sam.sort.bam
do
bamCoverage -p 10 -bs 20 -b ${file} -o ${file}.bw
done
sh 06.bam_to_bigwig.sh # bw文件可以在IGV中可视化
deeptools介绍：BAM文件转换为bigwig文件，转换为bigwig之前需要先对bam文件排序
Picard去重
# Picard下载地址：https://broadinstitute.github.io/picard/
vim 07.picard.sh
for file in *de-rRNA.fq.gz.sam.sort.bam
java -Xms20g -Xmx20g -XX:ParallelGCThreads=4 -jar
picard.jar MarkDuplicates I=${file} O=${file}.bam
M={file}.matrix ASO=coordinate
REMOVE_DUPLICATES=true 2> ${file}.picard.log &
sh picard.sh
peak鉴定
# peak calling
vim 08.macs2.sh
macs2 callpeak -t xx-1-ip_R2_trimmed.fq.gz_de-rRNA.fq.gz.sam.sort.bam xx-2-ip_R2_trimmed.fq.gz_de-rRNA.fq.gz.sam.sort.bam -c xx-1-input_R2_trimmed.fq.gz_de-rRNA.fq.gz.sam.sort.bam xx-2-input_R2_trimmed.fq.gz_de-rRNA.fq.gz.sam.sort.bam -f BAM -g 1.3e8 --nomodel -q 0.01 -n xx_peak
sh 08.macs2.sh
# 常用参数
-t/–TREATMENT FILENAME：treat组/IP 组
-c/–CONTROL：input组 或 mock(非特异性抗体，如IgG)组control:mock: 1）未使用抗体富集组；2）非特异性抗体，如IgG
-n/–NAME：为MACS2输出文件命名
-OUTDIR：MACS2结果文件保存路径
-f/–FORMAT FORMAT："SAM", "BAM", "BOWTIE", "BAMPE" or "BEDPE"；默认自动检测输入文件格式
-g/–GSIZE：有效基因组大小；基因组中有大量重复序列测序测不到，原基因组90% 或 70%；人类默认值是– 2.7e
-q/–QVALUE：qvalue 设定call significant regions的阈值；默认，0.01
Q-values are calculated from p-values using Benjamini-Hochberg procedure.
-p/–PVALUE：设定p值时， qvalue不再起作用。
-m/–-MFOLD：（DEFAULT:5, 50)
–SLOCAL, –LLOCAL：设定两个水平检测peak 区域，从而计算最大λ作为local λ。
默认，MACS2采用1000bp为small local region(–slocal)，10000bps为large local region(–llocal)区域设置的太小，
尖峰会掩盖掉旁边显著性的峰。
–nomodel：MACS2不构建模型。
--to-large 把小样本放大成和大样本一致
--down-sample. Macs2会对样品随机抽样，每次结果不一样（不建议使用）
--keep-dup 保留重复(基于课题设计)
macs2介绍：基于新的模型可以识别染色质免疫沉淀技术转录因子结合位点。macs 可以直接应用于ChIP-Seq 数据，也可以将ChIP-Seq数据与control结合起来提高特异性。
peak注释和motif预测
# Homer对peak进行注释http://homer.ucsd.edu/homer/ngs/annotation.html
# 下载配置文件
wget http://homer.salk.edu/homer/configureHomer.pl
# 使用配置文件进行软件配置
perl configureHomer.pl -install
# 下载hg19 基因组信息
perl configureHomer.pl -install hg19
vim 09.homer.sh
perl annotatePeaks.pl xx_peak hg19 -annStats ann_stats.log > xx_peak.ann.tsv
sh 09.homer.sh
head xx_peak.ann.tsv #注释结果文件
less ann_stats.log #peak分布统计情况
# motif finding
vim 10.motif.sh
perl findMotifsGenome.pl xx_peak hg19 xx_motif -len 8,10,12
sh 10.motif.sh
差异基因和差异peak
# 差异基因
vim 11.cutdiff.sh
cutdiff -o cutdiff_output -p 4 -L ctrl，treat.gtf xx-1-input_R2_trimmed.fq.gz_de-rRNA.fq.gz.sam.sort.bam xx-2-input_R2_trimmed.fq.gz_de-rRNA.fq.gz.sam.sort.bam
sh 11.cutdiff.sh
# 差异m6Apeak
vim exomepeak.sh
library ("exomepeak")
gtf=c(".gtf")
ctrl_IP=c("xx-1-ip_R2_trimmed.fq.gz_de-rRNA.fq.gz.sam.sort.bam")
ctrl_input=c("xx-1-input_R2_trimmed.fq.gz_de-rRNA.fq.gz.sam.sort.bam")
treat_IP=c("xx-2-ip_R2_trimmed.fq.gz_de-rRNA.fq.gz.sam.sort.bam")
treat_input=c("xx-2-input_R2_trimmed.fq.gz_de-rRNA.fq.gz.sam.sort.bam")
diff_peak <- exomepeak(GENE_ANNO_GTF=gtf,
IP_BAM=c(ctrl_IP),INPUT_BAM=c(ctrl_input),
TREATED_IP_BAM=c(treat_IP),TREATED_INPUT_BAM=c(treat_input),
EXPERIMENT_NAME=c("xx")
