散点图加密度图

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

基因表达趋势分析(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))

limma去除批次效应

什么是批次效应(batch effect)?
不同平台的数据,同一平台的不同时期的数据,同一个样品不同试剂的数据,以及同一个样品不同时间的数据等等都会产生一种batch effect 。这种影响如果广泛存在应该被足够重视,否则会导致整个实验和最终的结论失败。比对实验组和对照组,不同的处理是患病和不患病(测序时,先测得疾病,然后测得正常),然后你通过分析,得到很多差异表达的基因。现在问题来了,这个差异表达的结果是和你要研究的因素有关,还是时间有关,这个问题里时间就会成为干扰实验结果的因素,这个效应就是batch effect。

library(limma)
data <- read.table("mouse_gc.xls",sep="\t",header = TRUE)
batch1 <- c(rep('fish',9),rep('mouse',9))
batch1 <- as.factor(batch1)
design <- model.matrix(~0 + batch1)
#data=normalizeBetweenArrays(data)#如果组内的中位数不在同一条水平线上,现在该参数校正,再使用下面的批次校正 new_data <- removeBatchEffect(data[,2:19], batch = batch1) boxplot(data[,2:19]) boxplot(new_data) write.table(new_data, file="removeBatch.xls", sep="\t", quote=F, row.names=T, col.names=T)

limma结果出现负值请参考:https://cloud.tencent.com/developer/article/1680305
参考链接:https://www.omicsclass.com/article/1113

preprocessCore分位数标准化

library("preprocessCore")
library("ggplot2")
data <- read.table("mouse_gc.xls",sep="\t",header = TRUE)
rownames(data) <- data$GeneID
#rownames(data) <- data[,1]
#data_mat <- data.matrix(data[,-1]) 
quantile_normalisation <- function(df){
  df_rank <- apply(df,2,rank,ties.method="min")
  df_sorted <- data.frame(apply(df, 2, sort))
  df_mean <- apply(df_sorted, 1, mean)

  index_to_mean <- function(my_index, my_mean){
    return(my_mean[my_index])
  }

  df_final <- apply(df_rank, 2, index_to_mean, my_mean=df_mean)
  rownames(df_final) <- rownames(df)
  return(df_final)
}
new_data <- quantile_normalisation(data[,2:19])
boxplot(data)
boxplot(new_data)
write.table(new_data, file="111.xls", sep="\t", quote=F, row.names=T, col.names=T)
##############################################################
col_names <- colnames(data)
row_names <- rownames(data)
data <- as.matrix(data)
b=normalize.quantiles(data)
boxplot(b)
##使用normalize.quantiles函数和上面结果相同

数据相关性分析

简介:评价两组数据之间的相关性,有皮尔森(pearson)相关系数,斯皮尔曼(spearman)相关系数和肯德尔(kendall)相关系数。在这三大相关系数中,spearman和kendall属于等级相关系数亦称为“秩相关系数”,是反映等级相关程度的统计分析指标。相关系数的绝对值越大,相关性越强,相关系数越接近于1或-1,相关度越强,相关系数越接近于0,相关度越弱。

Pearson(皮尔逊)相关系数:皮尔逊相关也称为积差相关(或积矩相关)是英国统计学家皮尔逊于20世纪提出的一种计算直线相关的方法。适用于:

(1)、两个变量之间是线性关系,都是连续数据。

(2)、两个变量的总体是正态分布,或接近正态的单峰分布。

(3)、两个变量的观测值是成对的,每对观测值之间相互独立。

Spearman Rank(斯皮尔曼等级)相关系数:在统计学中,斯皮尔曼等级相关系数以Charles Spearman命名,并经常用希腊字母ρ(rho)表示其值。又称秩相关系数,是利用两变量的秩次大小作线性相关分析,对原始变量的分布不作要求,属于非参数统计方法,适用范围要广些。斯皮尔曼等级相关是根据等级资料研究两个变量间相关关系的方法。它是依据两列成对等级的各对等级数之差来进行计算的,所以又称为“等级差数法”。斯皮尔曼等级相关对数据条件的要求没有积差相关系数严格,只要两个变量的观测值是成对的等级评定资料,或者是由连续变量观测资料转化得到的等级资料,不论两个变量的总体分布形态、样本容量的大小如何,都可以用斯皮尔曼等级相关来进行研究。对于服从Pearson 相关系数的数据亦可计算Spearman 相关系数,但统计效能要低一些。Pearson 相关系数的计算公式可以完全套用 Spearman 相关系数计算公式,但公式中的x 和y 用相应的秩次代替即可。

Kendall Rank(肯德尔等级)相关系数:在统计学中,肯德尔相关系数是以Maurice Kendall命名的,并经常用希腊字母τ(tau)表示其值。用于反映分类变量相关性的指标,适用于两个分类变量均为有序分类的情况。对相关的有序变量进行非参数相关检验; 取值范围在-1-1 之间,此检验适合于正方形表格;肯德尔(Kendall)W 系数又称和谐系数,是表示多列等级变量相关程度的一种方法。适用这种方法的数据资料一般是采用等级评定的方法收集的,即让K 个 评委(被试)评定N 件事物,或1 个评委(被试)先后K 次评定N 件事物。等级评定法每个评价者对N 件事物排出一个等级顺序,最小的等级序数为1 ,最大的为N ,若并列等级时,则平分共同应该占据的等级,如,平时所说的两个并列第一名,他们应该占据1 ,2 名,所以它们的等级应是1.5, 又如一个第一 名,两个并列第二名,三个并列第三名,则它们对应的等级应该是1,2.5,2.5,5,5,5, 这里2.5 是2,3 的平均,5 是4,5,6 的平均。

R包计算相关性系数

library(ggplot2)
library(ggpubr)
data <- read.table("111.xls",sep="\t",header = TRUE)#列是特征(两个特征),行是物种名
ggplot(data=data, aes(x=genome, y=Tes))+geom_point(color="red")+stat_smooth(method="lm")+stat_cor(data=data, method = "pearson")+theme_classic()
#stat_cor(data=dat, method = "pearson")意为用pearson相关进行相关性分析,可以自行更改方法
#stat_smooth是画拟合曲线的函数
#se=FALSE意思为不画出置信区间
#data=后跟需要画图的数据的文件名
#X=后跟作为X轴的数据的那一列的列名
#Y=后跟作为Y轴的数据的那一列的列名
#geom_point函数是个性化设置散点图点的形状,颜色,大小等,此处只设置了颜色,有需要可自行加入
#theme_classic() x,y轴实线化
library(corrplot)
library(pheatmap)
corr <- cor(data, method = 'pearson')
pheatmap(corr)
corrplot(corr, type = 'upper', tl.col = 'black', order = 'hclust', tl.srt = 45, addCoef.col = 'white')