ViewBS:DNA甲基化数据的可视化

官网提供了三种安装方式:conda、docker和正常一步步安装,前两种较为简单,能解决较多的依赖关系,下边介绍一步步安装(https://github.com/xie186/ViewBS)

#安装htslib
git clone https://github.com/samtools/htslib.git
git submodule update --init --recursive
autoreconf
./configure
make
sudo make insatll
#安装ViewBS
wget -c https://github.com/xie186/ViewBS/archive/refs/tags/v0.1.10.tar.gz
tar -zxvf v0.1.10.tar.gz
cd ViewBS-0.1.10
/opt/biosoft/ViewBS-0.1.10/ext_tools/cpanm --local-lib=~/perl5 local::lib && eval $(perl -I ~/perl5/lib/perl5/ -Mlocal::lib)
/opt/biosoft/ViewBS-0.1.10/ext_tools/cpanm Getopt::Long::Subcommand
/opt/biosoft/ViewBS-0.1.10/ext_tools/cpanm Getopt::Long (> 2.50)
/opt/biosoft/ViewBS-0.1.10/ext_tools/cpanm --force Bio::DB::HTS::Tabix

ViewBS主要以Bismark的bismark_methylation_extractor输出结果为输入,同时根据不同的目的需要准备全基因组DNA序列的fasta文件,所有基因的bed文件,转座子bed文件,差异基因的bed文件等。

#数据准备
samtools faidx genome.fasta
bgzip ../A.1_bismark_bt2_pe.deduplicated.CX_report.txt ./ #Bismark结果需要bgzip压缩
tabix -p vcf A.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz #生成tbi结尾的index文件
#分析流程
ViewBS MethCoverage --reference /home/wuchangsong/gc_genome/17.Geta/0.initial_data/genome.fasta --sample A.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-1 --sample B.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-2 --sample C.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-3 --sample D.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-1 --sample E.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-2 --sample F.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-3 --sample G.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-1 --sample H.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-2 --sample I.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-3 --outdir MethCoverage --prefix BS_seq_allsam
ViewBS MethGeno --genomeLength /home/wuchangsong/gc_genome/17.Geta/0.initial_data/genome.fasta.fai --sample A.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-1 --sample B.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-2 --sample C.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-3 --sample D.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-1 --sample E.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-2 --sample F.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-3 --sample G.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-1 --sample H.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-2 --sample I.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-3 --prefix BS_geno_sample --context CG --outdir MethGeno_100k --minLength 1000000 --win 100000 --step 100000
ViewBS MethOverRegion --region repeats.bed --sample A.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-1 --sample B.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-2 --sample C.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-3 --sample D.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-1 --sample E.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-2 --sample F.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-3 --sample G.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-1 --sample H.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-2 --sample I.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-3 --prefix bis_gene_all_sample --context CG --outdir MethOverRegion

详细流程参考官网

散点图加密度图

library(ggplot2)
library(tidyverse)
library(ggExtra)
data = read.table("mulvscon_ma.xls",header=T,row.names=1)
aes = aes(x=log2(baseMean),y=log2FoldChange)
maplot = ggplot(data=data,aes) + geom_point(aes(color=significance),size=0.5, show.legend = F)
p1 = maplot + scale_color_manual(values=c("#a92323","#999999","#a92323")) + geom_hline(yintercept=0,linetype=2,color="black")+geom_hline(yintercept=0.5,linetype=4,color="black")+ geom_hline(yintercept=-0.5,linetype=4,color="black")+theme_classic()
p2=ggMarginal(p1, type="density", fill="purple", color="purple")
p2

Science:新研究揭示B1细胞的起源

在一项新的研究中,德国马克斯-德尔布吕克分子医学中心的Klaus Rajewsky教授及其团队报道B1细胞的产生并不需要不同的祖细胞。相反,他们的实验表明B1典型的B细胞受体(B1-typical B-cell receptor)能够将B2细胞重编程为B1细胞,这表明B1细胞是由于它们的特殊B细胞受体而出现的。这项研究可能会解决一项持续了几十年的免疫学争论。相关研究结果发表在2019年2月15日的Science期刊上,论文标题为“BCR-dependent lineage plasticity in mature B cells”。在与疾病的斗争中,有一种东西是至关重要的:B细胞。这些特殊的细胞属于一类称为淋巴细胞的白细胞,是免疫系统中唯一能够产生抗体的细胞。这些Y形蛋白(即抗体)结合到在细菌或病毒等外来结构上,从而将它们标记为入侵者,以便吞噬细胞和其他的免疫细胞清除它们。

B1细胞已存在于新生儿中并在天然免疫(也称先天免疫)中起重要作用

B细胞有两种类型。B2细胞构成体内白细胞群体的最大部分,主要在血液和胸腺、脾脏、淋巴结、骨髓等淋巴器官中循环。另一方面,B1细胞主要存在于腹膜腔和胸膜腔中,因此存在于肠道和肺部周围的区域。它们对广泛的外来蛋白(称为抗原)作出反应,但也对身体自身的一些抗原作出反应,因而不同于高度特化的B2细胞。

B1细胞占新生儿所有B淋巴细胞的大部分,但是在成年人中,B1细胞所占的比例仅下降到百分之几。这就是B1细胞被认为是天然免疫载体—先天免疫系统—而B2细胞主要负责适应性免疫的的原因之一。当遭受感染或接种疫苗后,适应性免疫就会出现。

几十年来,免疫学家们一直在争论B细胞的起源

在此之前,科学家们仍不清楚不同类型的B细胞是如何产生的。 论文共同第一作者、马克斯-德尔布吕克分子医学中心免疫调节与癌症研究小组成员Robin Graf博士说,“一些免疫学家认为B1细胞和B2细胞来自不同的祖细胞。其他免疫学家认为特殊的自身反应性B细胞受体会触发B1细胞形成。”

如今,这项新研究为第二个假设的有效性提供了明确的证据。Graf解释道,“我们利用仅在B1细胞中发现的B1典型B细胞受体替换成熟的B2细胞中的B细胞受体。”

经过操纵的B2细胞呈现出B1细胞的特性

这一过程将这些经过操纵的B2细胞转化为B1细胞。Graf报道,“我们能够证实这些细胞获得了B1典型的表面标志物。”这些经过操纵的B2细胞也具有B1细胞的功能特性。Graf 说,“当我们将它们移植到小鼠体内时,它们会归巢到身体中天然发现B1细胞的那些部位。”

此外,这些细胞开始自发地产生抗体。Graf解释道,“这也是B1细胞的典型特征。”更重要的是,一旦B1典型的B细胞受体在B2细胞表面上表达,这些细胞在一到两周的时间内开始大量增殖。这与B1细胞在早期阶段的自然发育非常相似—这个过程很少被研究过。

旷日持久的争论即将结束

在这项研究的后期,Graf测量了这些经过操纵的B2细胞中数千个基因的活性。Graf 报道,“在这项研究中,我们发现这些相同的基因在这些细胞中和正常的B1细胞中都是有活性的。”

Graf期待这项新的研究将让免疫学家相信B1细胞是由于它们的B细胞受体的特异性而出现的。(生物谷 Bioon.com)

转载自:http://news.bioon.com/article/6734027.html

基因表达趋势分析(TCseq使用)

TCseq可根据不同的聚类方法将基因按照表达模式分类

BiocManager::install("TCseq")
library(TCseq)
data <- read.delim('df_TCseq.xls', row.names = 1, sep = '\t', check.names = FALSE)
data <- as.matrix(data)
tca <- timeclust(data, algo = "cm", k = 6, standardize = TRUE)
#character string giving a clustering method. Options are km' (kmeans), 'pam' (partitioning around medoids), 'hc' (hierachical clustering), 'cm' (cmeans).
p <- timeclustplot(tca, value = "z-score(TPM)", cols = 3)#所有cluster合并一个图
print(p[[1]])#单一cluster作图
a <- as.data.frame(tca@cluster)
table(a)
names(a) <- 'Cluster'
Cluster2 <- subset(a,Cluster == 2)
write.table(Cluster2, file="Cluster2_gene.xls", sep="\t", quote=F, row.names=T, col.names=T)

导出不同cluster基因做功能富集

BS-seq分析流程(二)

DNA甲基化:DNA甲基化为DNA化学修饰的一种形式,能在不改变DNA序列的前提下,改变遗传表现。为表观遗传编码的一部分,是一种外遗传机制。DNA甲基化过程会使甲基添加到DNA分子上,例如在胞嘧啶环的5’碳上:这种5’方向的DNA甲基化方式可见于所有脊椎动物。 在人类细胞内,大约有1%的DNA碱基受到了甲基化。

gzip -dc A.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz > A.1_bismark_bt2_pe.deduplicated.CX_report.txt
perl BismarkCX2methykit.pl A.1_bismark_bt2_pe.deduplicated.CX_report.txt
#准备注释需要的bed文件(格式12)
convert2bed -i gtf < out.gtf > out6.bed
grep 'exon' out6.bed > aa && mv aa out6.bed
python3 bed6Tobed12.py out6.bed  > out12.bed#bed6Tobed12.py
#methylKit安装及使用
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("methylKit")
library(methylKit)
file.list = list("A.CG_methykit.txt","B.CG_methykit.txt","C.CG_methykit.txt","D.CG_methykit.txt","E.CG_methykit.txt","F.CG_methykit.txt","G.CG_methykit.txt","H.CG_methykit.txt","I.CG_methykit.txt")
myobjDB=methRead(file.list,
sample.id=list("Control_1","control_2","control_3","Single_1","Single_2","Single_3","Multiple_1","Multiple_2","Multiple_3"),
assembly="HZGC",
treatment=c(0,0,0,1,1,1,2,2,2),
context="CpG",
mincov = 10,
dbtype = "tabix",
dbdir = "methylDB",
)
filtered.myobj=filterByCoverage(myobjDB,lo.count=10,lo.perc=NULL, hi.count=NULL,hi.perc=99.9)
meth=unite(filtered.myobj, destrand=FALSE)

pdf("Correlation.pdf")
getCorrelation(meth,plot=TRUE)
dev.off()
pdf("clusterSamples.pdf")
clusterSamples(meth, dist="correlation", method="ward", plot=TRUE)
dev.off()
pdf("PCAscreeplot.pdf")
PCASamples(meth, screeplot=TRUE)
dev.off()
pdf("PCAcluster.pdf")
PCASamples(meth)
dev.off()

tiles=tileMethylCounts(filtered.myobj,win.size=1000,step.size=1000)
meth=unite(tiles, destrand=FALSE)
myDiff=calculateDiffMeth(meth,num.cores=10)
diffMethPerChr(myDiff,plot=FALSE,qvalue.cutoff=0.01, meth.cutoff=25)
myDiff25p.hyper=getMethylDiff(myDiff,difference=25,qvalue=0.01,type="hyper")
myDiff25p.hypo=getMethylDiff(myDiff,difference=25,qvalue=0.01,type="hypo")
myDiff25p=getMethylDiff(myDiff,difference=25,qvalue=0.01)
diffMethPerChr(myDiff,plot=FALSE,qvalue.cutoff=0.01, meth.cutoff=25)
save.image("MS_M_S.RData")
Diffdataframe=getData(myDiff25p)
write.table(Diffdataframe, file="Diffdataframe.xls", sep="\t", quote=F, row.names=T, col.names=T)

R根据基因的表达量筛选基因

data <- read.delim('mouse_gc_log2.xls', row.names = 1, sep = '\t', check.names = FALSE)#行为基因ID列为样本的表达矩阵
data[order(rowSums(data),decreasing=T)[1:3000],]  #筛选表达量高的前3000个基因
data[rowSums(data)>1,]  #过滤掉表达量低的基因

data1=rowMeans(data)>1  #按平均数筛选
data2=rowSums(data>0)>6 #表达量不为0的样品个数筛选
data=data[data1 & data2,] #联合一下
data$cv <- apply(data, 1, function(x){
  sd(x)/mean(x)*100
})
data_df <- data[order(data$cv, decreasing = T)[1:3000], 1:15]#筛选变异系数最大的3000个基因

Gene-rank score计算

data <- read.delim('mouse_gc.xls', row.names = 1, sep = '\t', check.names = FALSE)
col_names <- colnames(data)
row_names <- rownames(data)
data <- as.matrix(data)
b <- apply(data, 2, function(y) rank(y) / length(y))
write.table(b, file="gene_rank_score.xls", sep="\t", quote=F, row.names=T, col.names=T)
#gene-rank没有改变样本间的相关性

MetabolAnalyze实现Pareto scaling

BiocManager::install("MetabolAnalyze")
library("MetabolAnalyze")
data <- read.delim('mouse_gc_TMM_cross.xls', row.names = 1, sep = '\t', check.names = FALSE)
col_names <- colnames(data)
row_names <- rownames(data)
data <- as.matrix(data)
b = scaling(data,type = "pareto")
boxplot(b)
write.table(b, file="Pareto_scaled.xls", sep="\t", quote=F, row.names=T, col.names=T)

edgeR之TMM标准化

library(edgeR)
library(ggplot2)
data <- read.delim('../mouse_TPM_notCross.matrix.xls', row.names = 1, sep = '\t', check.names = FALSE)
group <- factor(rep(c('A', 'B', 'C'), each = 3))
y <- DGEList(counts = data, group = group)
y <- calcNormFactors(y)
logcpm <- cpm(y, prior.count=3, log=TRUE)
write.table(logcpm, file="mouse_TMM.xls", sep="\t", quote=F, row.names=T, col.names=T)
#dgList <- estimateCommonDisp(y)
#dgList <- estimateTagwiseDisp(dgList)
#norm_counts.table <- t(t(dgList$pseudo.counts)*(dgList$samples$norm.factors))
#write.table(norm_counts.table, file="mouse_gc_pareto_TMM.xls", sep="\t", quote=F)
#pheatmap(log(data+1),cluster_rows=T,cluster_cols=T,scale="none",border_color="white",color=colorRampPalette(rev(c("red","white","blue")))(102))