Tag Archives: 绘图

七 21

R语言用hclust进行聚类分析

Posted on by

聚类的基础就是算出所有元素两两间的距离，我们首先做一些示例数据，如下：

x=runif(10)

y=runif(10)

S=cbind(x,y)                                 #得到2维的数组

rownames(S)=paste("Name",1:10,"")             #赋予名称，便于识别分类

out.dist=dist(S,method="euclidean")           #数值变距离

这个代码运行得到的S是一个矩阵，如下

> S

x         y

Name 1   0.41517985 0.4697017

Name 2   0.35653781 0.1132367

Name 3   0.52253349 0.3680286

Name 4   0.80558684 0.9834687

Name 5   0.04564145 0.8560690

Name 6   0.11044397 0.2988598

Name 7   0.34984447 0.8515141

Name 8   0.28097709 0.1260050

Name 9   0.81771888 0.5976135

Name 10 0.40700158 0.5236567

可以看出里面共有10个点，它们的X,Y坐标均已知，我们有6总方法可以求矩阵

image001

注释：在聚类中求两点的距离有：

1，绝对距离：manhattan

2，欧氏距离：euclidean 默认

3，闵科夫斯基距离：minkowski

4，切比雪夫距离：chebyshev

5，马氏距离：mahalanobis

6，蓝氏距离：canberra

用默认的算法求出距离如下

算出距离后就可以进行聚类啦！

out.hclust=hclust(out.dist,method="complete") #根据距离聚类

注释：聚类也有多种方法：

1，类平均法：average

2，重心法：centroid

3，中间距离法:median

4，最长距离法：complete 默认

5，最短距离法：single

6，离差平方和法：ward

7，密度估计法：density

接下来把聚类的结果图画出来

plclust(out.hclust)                           #对结果画图

image003

rect.hclust(out.hclust,k=3)                   #用矩形画出分为3类的区域

image005

out.id=cutree(out.hclust,k=3)                 #得到分为3类的数值

这里的out.id就是把每个点都分类了的分类数组，1,2,3.

 

七 21

用R语言批量做T检验

Posted on by

需要做T检验的的数据如下：其中加粗加黑的是case，其余的是control，需要进行检验的是

p_ma    p_mg    p_ag    p_mag 这四列数据，在case和control里面的差异，当然用肉眼很容易看出，control要比case高很多

 

individual m a g ma mg ag mag p_ma p_mg p_ag p_mag
HBV10-1 33146 3606 2208 308 111 97 39 0.79% 0.29% 0.25% 0.10%
HBV15-1 23580 10300 3140 1469 598 560 323 4.19% 1.71% 1.60% 0.92%
HBV3-1 107856 26445 15077 4773 2383 1869 1130 3.34% 1.67% 1.31% 0.79%
HBV4-1 32763 6448 4384 579 291 295 124 1.35% 0.68% 0.69% 0.29%
HBV6-1 122396 6108 7953 911 796 347 144 0.67% 0.59% 0.26% 0.11%
IGA-1 31337 22167 14195 3851 2752 4101 2028 6.50% 4.65% 6.92% 3.42%
IGA-10 6713 9088 12801 2275 2470 4284 1977 10.54% 11.44% 19.85% 9.16%
IGA-11 61574 23622 15076 5978 4319 3908 2618 6.71% 4.84% 4.38% 2.94%
IGA-12 38510 23353 20148 6941 6201 6510 4328 10.38% 9.27% 9.73% 6.47%
IgA-13 155379 81980 65315 26055 20085 17070 10043 10.38% 8.00% 6.80% 4.00%
IgA-14 345430 86462 71099 27541 21254 10923 6435 6.06% 4.67% 2.40% 1.42%
IgA-15 3864 3076 1942 378 207 389 145 4.66% 2.55% 4.80% 1.79%
IgA-16 3591 4106 2358 427 174 424 114 4.64% 1.89% 4.61% 1.24%
IgA-17 893 1442 799 68 28 78 18 2.27% 0.94% 2.61% 0.60%
IGA-2 23097 5083 5689 910 951 1173 549 2.89% 3.02% 3.72% 1.74%
IGA-3 14058 9364 8565 2130 1953 2931 1436 8.03% 7.36% 11.05% 5.41%
IGA-4 81571 36894 33664 11346 10131 9908 6851 8.86% 7.91% 7.74% 5.35%
IGA-5 27626 6492 4503 963 752 942 410 2.64% 2.06% 2.58% 1.12%
IGA-7 55348 5956 4028 833 476 468 207 1.30% 0.74% 0.73% 0.32%
IGA-8 31671 17097 10443 3894 2666 3514 2003 7.56% 5.17% 6.82% 3.89%

 

但是我们需要写程序对这些百分比在case和control里面进行T检验。

a=read.table("venn_dat.txt",header=T)

type=c(rep("control",5),rep("case",12),"control","control","case")

for (i in 9:12)

{

dat=as.numeric(unlist(lapply(a[i],function(x) strsplit(as.character(x),"%"))))

filename=paste(names(a)[i],".png",sep="")

png(file=filename, width = 320, height = 360)

b=t.test(dat~type)

txt=paste("p-value=",round(b$p.value[1],digits=4),sep="")

plot(as.factor(type),dat,ylab="percent(%)",main=names(a)[i],sub=txt,cex.axis=1.5,cex.sub=2,cex.main=2,cex.lab=1.5)

dev.off()

}

得到的图片如下：

image008

七 21

用R语言批量画韦恩图

Posted on by

需要画韦恩图的文件如下所示：

#CDR3_aa    count_all    count_IgM    count_IgA    count_IgG
CARGVDAGVDYW    243    25    196    22
CARHPRNYGNFDYW    174    171    3    0
CARENTMVRGVINPLDYW    166    8    75    83
CAREASDSISNWDDWYFDLW    129    15    114    0
CARDPDNSGAFDPW    118    1    117    0
CAKDLGGYW    98    3    4    91
CAREVADYDTYGWFLDLW    95    26    68    1
CVRNRGFFGLDIW    82    0    1    81
CARRSTNYHGWDYW    80    3    2    74

此处省略一万行。

简单解释一下数据，第一列是CDR3序列，我们需要对count_IgM    count_IgA    count_IgG这三列数据进行画韦恩图，数字大于0代表有，数字为0代表无。

这样我们根据序列就能得出每列数据所有的CDR3序列，即不为0的CDR3序列

每个个体都会输出一个统计文件，共20个文件需要画韦恩图

image005

对这个统计文件就可以进行画韦恩图。

画韦恩图的R代码如下：

library(VennDiagram)

files=list.files(path = ".", pattern = "type")

for (i in files){

a=read.table(i)

individual=strsplit(i,"\\.")[[1]][1]

image_name=paste(individual,".tiff",sep="")

IGM=which(a[,3]>0)

IGA=which(a[,4]>0)

IGG=which(a[,5]>0)

venn.diagram(list(IGM=IGM,IGA=IGA,IGG=IGG), fill=c("red","green","blue"), alpha=c(0.5,0.5,0.5), cex=2, cat.fontface=4, fontfamily=3, filename=image_name)

}

image007

但事实上，这个韦恩图很难表现出什么，因为我们的每个个体的count_IgM    count_IgA    count_IgG总数不一样。