Category: R

twilight处理pangenome结果

github下载脚本（https://github.com/ghoresh11/twilight）需要的R包data.table，ggplot2，optparse Rscript classify_genes.R -p gene_presence_absence.Rtab -g groups.tab 必须参数： -p, –presence_absence：Roary 或 Panaroo输出的gene_presence_absence.Rtab文件 -g, –grouping：一个制表符分隔的分组文件可选参数： -o, –out：输出目录名称（默认=“out”）。 -m, –min_size：忽略少于min_size基因组的组（默认 = 10）。 -c, –core_threshold：用于定义每个组内的核心基因的阈值（默认值 = 0.95）。 -r, –rare_threshold：用于定义每组内稀有基因的阈值（默认值 = 0.15）。 -h, –help：打印帮助信息并退出。输出： 1、classification.tab 基因簇在每个分类中存在的形式（核心、辅助、和稀有） 2、frequencies.csv 基因簇在每个分离中的频率 3、plots 文件夹，包含柱状图和PCA图

2022年9月19日
R添加点的标签

geom_text_repel()和geom_label_repel() 参数： size = 4, #注释文本的字体大小 box.padding = 0.5, #字到点的距离 point.padding = 0.8, #字到点的距离，点周围的空白宽度 min.segment.length = 0.5, #短线段可以省略 segment.color = “black”, #segment.colour = NA, 不显示线段 nudge_x/y:数据点与相应数据标签的距离 segment.size:指定线段的粗细 arrow：angle: 箭头的尖角的角度；length: 箭头尖角的长度；ends: “last”, “first”, “both”, 指定线段的那端画箭头；type: “open”和”closed” 指定箭头是否为封闭的三角形。 force：重叠文本标签之间的排斥力，默认为 1。 force_pull：文本标签与其对应数据点之间的吸引力，默认为 1。 max.time：尝试解决重叠的最大秒数。默认为 0.5。 max.iter：尝试解决重叠的最大迭代次数。默认为 10000。 max.overlaps：排除重叠太多内容的文本标签。默认为 10。 xlim, ylim：x 和 y 轴的限制。文本标签将受到这些限制。…

2022年1月17日
散点图加密度图

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

2021年5月30日
基因表达趋势分析（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)”,…

2021年4月14日
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”,…

2021年4月8日
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个基因

2021年3月30日
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没有改变样本间的相关性

2021年3月22日
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)

2021年3月20日
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”,…

2021年3月18日
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

2021年3月16日