0.软件安装略
1.比对得到bam文件
bwa index ref.fa #个人感觉bwa比bowtie2在短reads的比对上更好
bwa mem -t 10 -R '@RG\tID:sp1\tPL:illumina\tLB:sp1\tSM:sp1' ref.fa 1.fq.gz 2.fq.gz | samtools view --threads 20 -Sb - | samtools sort -@ 36 -m 2G -O bam - -o sp1.sort.bam #对bam文件进行排序
#sp1改成你相应的样本名即可
去重复方法二选一即可
sambamba markdup -r sp1.sort.bam sp1.rmdup.bam #去除PCR重复， -r是去除重复序列，不加就是标记
#这个简单，用这个，而且这个还快，用这个
gatk MarkDuplicates -I sp1.sort.bam -O sp1.rmdup.bam -M tmp1 #去除PCR重复
#gatk用的picard，得加这个参数 --REMOVE_SEQUENCING_DUPLICATES true
#不然就只是标记重复序列但是不删除, 不想看命令行了
#明显不如sambmba简单，果断弃；实际证明运行速度也拉跨
samtools view -q 30 -f 2 -F 4 -O bam -o sp1.final.bam sp1.rmdup.bam #保留唯一比对，过滤低质量比对
samtools index sp1.final.bam #记得建索引
#如果有多个实验处理，比对到同一个参考基因组即可
- ATAC signal plot
awk '$3 == "gene"' ref.gff | awk 'BEGIN{FS="\t|=|;";OFS="\t"}{print $1,$4-1,$4}'> ref.gene.TSS.bed
#想看ATAC signal在转录起始位点的分布, 所以提取TSS的bed
#实际上就是注释到的gene的第一个碱基
#这里可以自定义为任何bed文件
bamCoverage --bam sp1.final.bam -o sp1.bw --effectiveGenomeSize 2471489134 #--normalizeUsing RPKM 不建议标准化
#--effectiveGenomeSize的获取参考https://zhuanlan.zhihu.com/p/437993675
#faCount ref.fa > ref.count #统计所有的N
#本质是genome length减去所有的N
computeMatrix reference-point --referencePoint TSS -p 15 -b 3000 -a 3000 -R ref.gene.TSS.bed -S sp1.bw \
--skipZeros --missingDataAsZero --outFileSortedRegions regions.genes.bed -o sp1.matrix.TSS.gz \
#https://www.jianshu.com/p/dc0c49bfcd41 参数参考这里
#最好自己用-h把所有参数都看一遍，这很重要
#这个 .bw
plotHeatmap -m sp1.matrix.TSS.gz -out sp1.ATAC.signal.pdf --plotFileFormat pdf #图+1
2.1.这个图是用gatk标记重复序列但是没有去除，并且没有保留唯一比对的结果
为什么TSS这里不是最高峰？？
2.2.用samtools去除重复序列，且保留唯一比对的结果
2.3. 这种图可以衍生画无数个图
eg: 基因数量在TAD边界上或者随便什么的分布
需要输入两个bed文件，第一个是基因数量的bed文件，将bed文件转为bw格式需要两步:
#bedtools genomecov -i input.bed -g chrom.sizes.bed -bg > XX.bedgraph
#这一步存疑
#XX.bedgraph 有四列, 染色体/起始bp/终止bp/value
#所以直接到这一步也可以似乎
bedGraphToBigWig XX.bedgraph chrom.sizes XX.bw
这里的XX.bw其实本质上就是ATAC signal的格式
第二个是TAD边界的bed文件或者是别的什么随便的bed文件：site.bed
同上述流程中的ref.gene.TSS.bed
computeMatrix reference-point --referencePoint TSS -p 15 -b 3000 -a 3000 -R site.bed -S XX.bw -o XX.matrix.site.gz --skipZeros --missingDataAsZero --outFileSortedRegions XX.bed #参数的解释这里不再赘述
plotHeatmap -m XX.matrix.site.gz -out out.pdf --plotFileFormat pdf
3.Call peaks
macs2 callpeak -t sp1.final.bam -f BAMPE -nomodel -g 2471489134 -n sp1.peak -B -q 0.0
网友评论