日常瞎掰
这一波疫情的反弹给上海很多地方按下了“暂停键”,从干饭人
变成 “饭人” ,也许只有一个通告的距离。一觉醒来,也许小区大门口就贴上了“封闭管理的告示”。本以为就像告示说的那样,“封闭48小时,做两次核酸检测”,就完事了。可事情并原本预设的那样顺利,一轮核酸检测结束,接着第二轮检测,第二轮还没有结束,更糟糕的事情又来了,本来只是出不了小区,但不知啥时候冒出一个“密接”让我们出不了楼道的大门了。这一下活生生一个“饭人”样子了,连迈出屋门的自由都失去了。坚守在家的这些时间才体会到古人说的一句话,“生命诚可贵,爱情价更高,若为自由故,二者皆可抛”有了更深的认识,真的“自由”确实是可望而不及的!“大敌当前”,我们仍需咬紧牙关齐心协力、众志成城,希望早日散去这满城的阴霾。
HiCExplorer
这些时日一直在做HiC-seq数据分析的事,当然这个技术也不是特别新鲜的事物,所以也有不少的工具可以用来分析,今天我们就来简单介绍一款——HiCExplorer。该工具是基于python编写,功能还是比较全面,从原始fastq可以得到交互矩阵,进而做TAD、Loop、A/B compartment这些分析内容。其他软件本人也没有接触过,但使用过这个软件后还是觉得挺好用的。下面一起来看看。
1、build matrix
双端的测序数据比对时,Read1和Read2分开比对,然后用于构建交互矩阵:
# mapping
bwa mem -A1 -B4 -E50 -L0 index_path \
mate_R1.fastq.gz 2>>mate_R1.log | samtools view -Shb - > mate_R1.bam
bwa mem -A1 -B4 -E50 -L0 index_path \
mate_R2.fastq.gz 2>>mate_R2.log | samtools view -Shb - > mate_R2.bam
# build
hicBuildMatrix --samFiles mate_R1.bam mate_R2.bam \
--binSize 10000 \
--restrictionSequence GATC \
--danglingSequence GATC \
--restrictionCutFile cut_sites.bed \
--threads 4 \
--inputBufferSize 100000 \
--outBam hic.bam \
-o hic_matrix.h5 \
--QCfolder ./hicQC
# Correction
hicCorrectMatrix diagnostic_plot -m hic_matrix.h5 -o hic_corrected.png
# Visualization
hicPlotMatrix -m hic_corrected.h5 -o hic_plot.png --log1p
交互矩阵展示如下:
在展示所有染色体的交互热图时,最好使用分辨率值大一点的,如1M的分辨率,因为bin太小的话,整个矩阵的交互关系对就会很多,需要消耗很大的资源,这时可使用hicMergeMatrixBins
命令可以合并矩阵的bin,从而达到节省运行资源的目的。
2、TAD calling
TAD结果展示时,内核还是用的pyGenomeTracks
软件的功能,所用的配置文件track.ini
里面的具体参数这里就不介绍了:
hicFindTADs -m hic_corrected.h5 --outPrefix hic_corrected --numberOfProcessors 16
hicPlotTADs --tracks track.ini --region chrX:99974316-101359967 --outFileName public.tabplot.pdf
展示如下:
3、A / B compartment
A/B compartment的计算过程并不是很复杂,最终得到PCA的PC1当作结果,此原始结果的正负并不能直接分别指代A、B成分,需要根据借助两种区域对应的基因表达量,组蛋白修饰等情况来矫正,然后定义A、B区域。A compartment, 对应开放染色质区域;而负数对应区域包含的基因个数较少,含量也低,将其定义为B compartment, 对应封闭染色质区域。示意图如下:
hicPCA -m hic_corrected.h5 -noe 2 -pm hic.cor.mat.h5 -f bigwig -o hic.pc1.bw hic.pc2.bw
hicPlotMatrix -m public.cor.mat.h5 –bigwig public.pc1.bw -o public.pca1.pdf
展示如下:
小问题提示
HiCExplorer虽然挺好用的,但python版本的软件(version2.2.1.1)还是有些小问题的,比如说hicDetectLoops
子命令就无法正常使用,其实作者已经在软件的模块中实现了这个功能,当前版本的所有模块都在软件包的安装目录下可以找到,如下所示:
既然功能都实现了,为什么还无法使用呢?这个当然应该是软件编写时的疏忽了,导致这个模块没有被加入到子命令里面,从而不可用。如果想手动添加,可到python的bin目录添加相应的子命令,具体可以参考已经存在的相关子命令。
结束语
作者能写出这样好用且功能挺全面的软件,还是挺厉害的,不过还是希望软件作者能够意识到,可以尽快修复这样的小BUG
,让大家可以用的更顺畅。如果大家想更详细的了解可以移步官方介绍: HiCExplorer。