一
读取数据
1.读取绝对丰度数据In [1]:df head(df,10)
Out[1]:
C-1 C-2 C-3 C-4 C-5 C-6 C-7 M-1 M-2 M-3 ... M-5 M-6 M-7 M-8 Abundance superkingdom phylum class order familyNo_Rank 32729 43885 36532 29335 21815 26879 41869 31030 29862 33787 ... 43164 30980 21853 21792 471140 - - - - -Lactobacillus 24438 24449 3760 18204 28795 4537 12831 24450 17354 26211 ... 30850 19648 21548 18263 291188 Bacteria Firmicutes Bacilli Lactobacillales LactobacillaceaeUnassigned 25317 21725 13405 19347 14697 16676 19162 15819 22083 21072 ... 12032 7750 33000 26360 279560 Bacteria Bacteroidetes Bacteroidia Bacteroidales PorphyromonadaceaeAlistipes 5720 6196 18812 14525 19759 21996 10696 4411 5670 3663 ... 5744 6580 5415 5571 146280 Bacteria Bacteroidetes Bacteroidia Bacteroidales RikenellaceaeParabacteroides 6092 2683 8164 10246 10619 8742 7885 4493 12209 2423 ... 5628 4943 6382 11633 118324 Bacteria Bacteroidetes Bacteroidia Bacteroidales PorphyromonadaceaeBacteroides 9803 3257 7719 5919 8039 6959 3502 6514 8694 3840 ... 7839 4368 11986 15711 116935 Bacteria Bacteroidetes Bacteroidia Bacteroidales BacteroidaceaeBarnesiella 2982 1414 2341 9184 3707 6875 3812 3183 2773 5250 ... 3128 6584 3764 6680 64267 Bacteria Bacteroidetes Bacteroidia Bacteroidales PorphyromonadaceaeEisenbergiella 680 2283 4357 1743 1251 6988 3731 1832 2204 1598 ... 1193 649 3466 4671 45890 Bacteria Firmicutes Clostridia Clostridiales LachnospiraceaeOscillibacter 653 714 4704 1569 2210 4680 4303 2748 4912 800 ... 2424 729 1287 2486 37738 Bacteria Firmicutes Clostridia Clostridiales RuminococcaceaeTuricibacter 1373 2366 20 1405 56 6 2091 5060 9 9229 ... 1 1 4360 1693 27673 Bacteria Firmicutes Erysipelotrichia Erysipelotrichales Erysipelotrichaceae
2.读取分组数据In [2]:
info = read.table('sample_groups.xls',header = F,col.names = c('sample','group'))head(info)
Out[2]:
sample groupC-1 ControlC-2 ControlC-3 ControlC-4 ControlC-5 ControlC-6 Control
3.检查样品是否正确In [3]:
info$sample
Out[3]:
C-1 C-2 C-3 C-4 C-5 C-6 C-7 M-1 M-2 M-3 M-4 M-5 M-6 M-7 M-8
4.根据分组文件,从丰度数据中提取样品的丰度信息In [4]:
df = df[as.character(info$sample)]head(df)
Out[4]:
C-1 C-2 C-3 C-4 C-5 C-6 C-7 M-1 M-2 M-3 M-4 M-5 M-6 M-7 M-8No_Rank 32729 43885 36532 29335 21815 26879 41869 31030 29862 33787 25628 43164 30980 21853 21792Lactobacillus 24438 24449 3760 18204 28795 4537 12831 24450 17354 26211 15850 30850 19648 21548 18263Unassigned 25317 21725 13405 19347 14697 16676 19162 15819 22083 21072 11115 12032 7750 33000 26360Alistipes 5720 6196 18812 14525 19759 21996 10696 4411 5670 3663 11522 5744 6580 5415 5571Parabacteroides 6092 2683 8164 10246 10619 8742 7885 4493 12209 2423 16182 5628 4943 6382 11633Bacteroides 9803 3257 7719 5919 8039 6959 3502 6514 8694 3840 12785 7839 4368 11986 15711
5.检查分组信息In [5]:
info$group
Out[5]:
Control Control Control Control Control Control Control NNK-BaP NNK-BaP NNK-BaP NNK-BaP NNK-BaP NNK-BaP NNK-BaP NNK-BaP
分组信息发现,仅有2个分组,组名分别是Control、NNK-BaP
二
计算该组的丰度之和、均值、方差、标准差
查看Control组包含哪些样品
In [6]:
info[info$group=='Control',1]
Out[6]:
C-1 C-2 C-3 C-4 C-5 C-6 C-7
apply函数计算Control组的和、均值、方差、标准差
In [7]:
df['Control_sum'] = apply(df[as.character(info[info$group=='Control',1])],1,sum)df['Control_mean'] = apply(df[as.character(info[info$group=='Control',1])],1,mean)df['Control_var'] = apply(df[as.character(info[info$group=='Control',1])],1,var)df['Control_sd'] = apply(df[as.character(info[info$group=='Control',1])],1,sd)head(df)
Out[7]:
C-1 C-2 C-3 C-4 C-5 C-6 C-7 M-1 M-2 M-3 M-4 M-5 M-6 M-7 M-8 Control_sum Control_mean Control_var Control_sdNo_Rank 32729 43885 36532 29335 21815 26879 41869 31030 29862 33787 25628 43164 30980 21853 21792 233044 33292.000 64182849 8011.420Lactobacillus 24438 24449 3760 18204 28795 4537 12831 24450 17354 26211 15850 30850 19648 21548 18263 117014 16716.286 99804027 9990.197Unassigned 25317 21725 13405 19347 14697 16676 19162 15819 22083 21072 11115 12032 7750 33000 26360 130329 18618.429 16946400 4116.601Alistipes 5720 6196 18812 14525 19759 21996 10696 4411 5670 3663 11522 5744 6580 5415 5571 97704 13957.714 43482964 6594.161Parabacteroides 6092 2683 8164 10246 10619 8742 7885 4493 12209 2423 16182 5628 4943 6382 11633 54431 7775.857 7342272 2709.663Bacteroides 9803 3257 7719 5919 8039 6959 3502 6514 8694 3840 12785 7839 4368 11986 15711 45198 6456.857 5800759 2408.477
apply函数计算NNK-BaP组的和、均值、方差、标准差
In [8]:
df['NNK-BaP_sum'] = apply(df[as.character(info[info$group=='NNK-BaP',1])],1,sum)df['NNK-BaP_mean'] = apply(df[as.character(info[info$group=='NNK-BaP',1])],1,mean)df['NNK-BaP_var'] = apply(df[as.character(info[info$group=='NNK-BaP',1])],1,var)df['NNK-BaP_sd'] = apply(df[as.character(info[info$group=='NNK-BaP',1])],1,sd)head(df)
Out[8]:
C-1 C-2 C-3 C-4 C-5 C-6 C-7 M-1 M-2 M-3 ... M-7 M-8 Control_sum Control_mean Control_var Control_sd NNK-BaP_sum NNK-BaP_mean NNK-BaP_var NNK-BaP_sdNo_Rank 32729 43885 36532 29335 21815 26879 41869 31030 29862 33787 ... 21853 21792 233044 33292.000 64182849 8011.420 238096 29762.000 48868388 6990.593Lactobacillus 24438 24449 3760 18204 28795 4537 12831 24450 17354 26211 ... 21548 18263 117014 16716.286 99804027 9990.197 174174 21771.750 25821419 5081.478Unassigned 25317 21725 13405 19347 14697 16676 19162 15819 22083 21072 ... 33000 26360 130329 18618.429 16946400 4116.601 149231 18653.875 72916680 8539.126Alistipes 5720 6196 18812 14525 19759 21996 10696 4411 5670 3663 ... 5415 5571 97704 13957.714 43482964 6594.161 48576 6072.000 5639229 2374.706Parabacteroides 6092 2683 8164 10246 10619 8742 7885 4493 12209 2423 ... 6382 11633 54431 7775.857 7342272 2709.663 63893 7986.625 22692800 4763.696Bacteroides 9803 3257 7719 5919 8039 6959 3502 6514 8694 3840 ... 11986 15711 45198 6456.857 5800759 2408.477 71737 8967.125 17710582 4208.394
三
按列计算相对丰度
In [9]:
df_Relative = apply(df,2,function(d) d/sum(d))head(df_Relative)
Out[9]:
C-1 C-2 C-3 C-4 C-5 C-6 C-7 M-1 M-2 M-3 ... M-7 M-8 Control_sum Control_mean Control_var Control_sd NNK-BaP_sum NNK-BaP_mean NNK-BaP_var NNK-BaP_sdNo_Rank 0.13115102 0.17585513 0.14639033 0.11755065 0.08741665 0.10770901 0.16777666 0.12434282 0.11966244 0.135390620 ... 0.08756892 0.08732449 0.13340707 0.13340707 0.24318010 0.14420668 0.11926172 0.11926172 0.19457417 0.10579967Lactobacillus 0.09792749 0.09797157 0.01506700 0.07294672 0.11538677 0.01818058 0.05141614 0.09797557 0.06954062 0.105032218 ... 0.08634673 0.07318314 0.06698518 0.06698518 0.37814391 0.17982493 0.08724334 0.08724334 0.10281045 0.07690602Unassigned 0.10144980 0.08705600 0.05371626 0.07752693 0.05889354 0.06682375 0.07678560 0.06338959 0.08849058 0.084439315 ... 0.13223697 0.10562929 0.07460741 0.07460741 0.06420761 0.07409938 0.07474945 0.07474945 0.29032475 0.12923607Alistipes 0.02292107 0.02482849 0.07538309 0.05820430 0.07917789 0.08814195 0.04286081 0.01767567 0.02272072 0.014678304 ... 0.02169888 0.02232400 0.05593109 0.05593109 0.16475105 0.11869583 0.02433160 0.02433160 0.02245313 0.03594017Parabacteroides 0.02441175 0.01075127 0.03271462 0.04105758 0.04255225 0.03503078 0.03159662 0.01800426 0.04892367 0.009709399 ... 0.02557383 0.04661554 0.03115927 0.03115927 0.02781887 0.04877431 0.03200385 0.03200385 0.09035356 0.07209653Bacteroides 0.03928239 0.01305139 0.03093143 0.02371850 0.03221373 0.02788597 0.01403315 0.02610278 0.03483843 0.015387575 ... 0.04803007 0.06295682 0.02587379 0.02587379 0.02197829 0.04335292 0.03593289 0.03593289 0.07051638 0.06369227
四
写入文件进行保存
In [10]:
write.csv(df_Relative,'genus.taxon.RelativeAbundance.cvs')
往期相关链接:
1、R基础篇
如何使用Rstudio练习R基础教程;
R相关软件及R包安装;
【零基础学绘图】之绘制venn图(五);
【零基础学绘图】之绘制barplot柱状图图(四);
【零基础学绘图】之绘制heatmap图(三);
【零基础学绘图】之绘制PCA图(二);
【零基础学绘图】之alpha指数箱体图绘制(一);
2、R进阶
【绘图进阶】之交互式可删减分组和显示样品名的PCA 图(三);
【绘图进阶】之绘制PCA biplot图(二);
【进阶篇绘图】之带P值的箱体图、小提琴图绘制(一);
3、数据提交
3分钟学会微生物多样性云平台数据分析;
3分钟学会CHIP-seq类实验测序数据可视化 —IGV的使用手册;
10分钟搞定多样性数据提交,最快半天内获取登录号,史上最全的多样性原始数据提交教程;
20分钟搞定GEO上传,史上最简单、最详细的GEO数据上传攻略;
4、表达谱分析
如何对GEO数据进行差异分析;
5、医学数据分析
【WGS服务升级】人工智能软件SpliceAI助力解读罕见和未确诊疾病中的非编码突变;
隐性疾病trio家系别忽视单亲二倍体现象——天昊数据分析助力临床疾病诊断新添UPD(单亲二倍体)可视化分析工具;
【昊工具】Oh My God! 太好用了吧!疾病或表型的关键基因查询数据库,我不允许你不知道Phenolyzer;
如果您对本文案介绍的方法或代码有疑问,
请扫码添加QQ群沟通
【本群将为大家提供】
分享生信分析方案
提供数据素材及分析软件支持
定期开展生信分析线上讲座
QQ号:1040471849
作者:大熊
审核:有才
来源:天昊生信团
有问题,找生信团,生信团帮您出谋划策,通用需求,我们会出专门教程哦。