Matlab error when running GISTIC

If you instal MCR (MATLAB Compiler Runtime) provided by GISTIC package, may have the following error. This error could disrupt GISTIC.
libGL error: failed to load driver: swrast

If this situation occurs, rename the file found at $MATLAB_ROOT/sys/os/glnxa64/libstdc++.so.6″ to “libstdc++.so.6.old”, This forces MATLAB to use the OS library.

Works for me.

Ref:
https://ww2.mathworks.cn/matlabcentral/answers/296999-libgl-error-unable-to-load-driver-in-ubuntu-16-04-while-running-matlab-r2013b

GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers

Prepare a data frame for sample CNV data

If we want to cluster samples based on CNV data, a dataframe is needed. However, CNV segments in each sample are not the same. Maybe overlap or distinct. I think CNTools package migh solve this challenge. An example is shown as below. The result is a reduced segment data frame.

BiocManager::install("CNTools")
data("sampleData")
seg <- CNSeg(sampleData)
rdseg <- getRS(seg, by = "region", imput = FALSE, XY = FALSE, what = "mean") 
View(rdseg@rs)

Input dataframe has six columns (“ID”,”chrom”,”loc.start”,”loc.end”,”num.mark”,”seg.mean”) including 277 samples and 54825 segments.

The result can be got from rdseg@rs, like this

Cheers

Also, we can use CNRegions from iClusterPlus package.

CNregions(sampleData)

Ref: https://www.rdocumentation.org/packages/CNTools

https://rdrr.io/bioc/iClusterPlus/man/CNregions.html

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

对Autoencoder(自编码器)的理解

通常数据的维度太大,可视化很难,也不利用模型的学习。有时候拿到数据做个PCA或者tSNE,就是把维度缩小到2维(当然也可以3维),便于看数据之间的关系。在机器学习中,Autoencoder也是一种降维的方式, Autoencoder输入层的神经元的数目和输出层的神经元的数目必须,而且要保证输出的结果尽最大可能和输入的结果一致。

图片来自网络

如上图所示,维度由大到小是decode过程,输出的结果可以从中间层经过encode得到,那么中间层保留了输入层的信息(因为输出层的结果从中间层得到),那么中间层的数据结果,就是降维后的结果,可以拿来做其他事情。 网络的复杂程度根据样本数设计。

无监督的聚类,便可以从中间层开始;数据的学习也可以从中间层开始。当输入层是多组学数据时,中间层便是融合后的结果。

FPKM转TPM

R code

fpkm2tpm = function(fpkm){
  exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}
tpm = apply(expMatrix, 2, fpkm2tpm)

If the expression matrix has NA value

fpkm2tpm <- function(fpkm){
  tpm <- exp(log(fpkm) - log(sum(fpkm,na.rm=T)) + log(1e6))
  tpm[which(is.na(tpm))] <- 0
  return(tpm)
}

TPMi=( FPKMi / sum(FPKMj ) * 10^6

可变多聚腺苷酸化Alternative Polyadenylation (APA) 检测

可变多聚腺苷酸化Alternative Polyadenylation (APA),如下图所示(图片来自参考),在不同的APA信号位点切割,然后添加polyA。这种调控机制属于转录后调控,可能会影响蛋白的序列(发生在编码区),也可能影响蛋白的稳定性(比如非编码区内的miRNA的调控区域)。其实也是可变剪接的一种情况。

常用的软件是Dapars,这个软件现在也有了升级的版本Dapars2。参考:

https://github.com/ZhengXia/dapars

https://github.com/3UTR/DaPars2

分析流程很相似,Dapars2多了 normalize library sizes 。

第一步:生成Wiggle文件

这一步其实是统计每个点或者区域的覆盖度,便于后续计算,用bedtools统计STAR比对后的结果即可。

genomeCoverageBed -bg -ibam Aligned.sortedByCoord.out.bam -split >  Aligned.sortedByCoord.out.wig

第二步: Generate region annotation: DaPars_Extract_Anno.py

这一步需要从bed12格式的文件中,提取远端ployA的位点,进而便于推测近端polyA的位点。

1)作者推荐从UCSC上下载,http://genome.ucsc.edu/cgi-bin/hgTables?command=start,选择bed格式之后,选这whole gene即可。

2)如果你有特定的gtf文件,可以将gtf文件转成bed12格式,先将gtf转成genPred格式,然后再转成bed12格式。http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/上有现成的工具

第三步(Dapars2): Generate mapped reads files for all Samples

这一步是为了 normalize library sizes ,统计每个样本比对上的reads数据,个是如下,与下一步的顺序一致。Dapars2才有这一步。

第四步: Run DaPars

配置文件,运行即可。结果中会告诉ΔDPUI也就是远端APA位点的利用率的差值及显著性,根据这些挑选候选。然后回到IGV上进行查看,如下基因方向从右到左,可以看出来前六个和后三个样本在最后一个外显子上的覆盖度截然不同,不同的地方可能就是发生了APA事件,然后添加polyA。

参考

Weil T T. Post-transcriptional regulation of early embryogenesis[J]. F1000prime reports, 2015, 7.

https://hpc.oit.uci.edu/~leil22/DaPars2_Documentation/DaPars2.html

https://www.jianshu.com/p/21b697cec428

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################