简介面向生物信息学与医学影像分析场景这套Matlab代码包用于将MRI脑图与来自微阵列数据的基因表达艾伦人类大脑图谱AHBA进行关联分析解决跨模态脑数据对齐与特征提取的问题。代码支持MATLAB 2014/2019a/2021a多种版本采用参数化编程风格关键参数可方便更改并附有案例数据支持直接运行适合计算机、电子信息工程、数学等专业的学生在课程设计、期末大作业和毕业设计中参考使用。包内共4个文件包含3个.m脚本和1个Markdown说明文档整个压缩包仅24KB结构非常精简。其中脚本分别承担数据读取与预处理、AHBA网络可视化、基于偏最小二乘的关联建模等任务代码注释明细、思路清晰README文档则对使用方法做了简要说明。目前已有119人学习下载对于需要快速掌握MRI与基因表达联合分析方法的同学和研究者这是一套可直接运行、易于二次开发的实用参考代码。1. 把MRI脑图对齐到艾伦图谱转录组-影像关联到底在算什么做神经影像的人手里常有一堆脑图皮层厚度、任务激活、静息态功能连接每一张都是全脑体素级别的空间模式。但光有模式不够大家真正想问的是“这种空间分布背后的分子基础是什么”。把MRI脑图与艾伦人类大脑图谱AHBA的微阵列基因表达数据关联起来就是回答这个问题的标准做法在MNI152空间里把两百多个脑区采样点的影像值取出来与两万多个基因的表达逐一算相关找出哪些基因的空间表达模式和你这张MRI图最像。这套MATLAB代码把坐标换算、探针聚合、相关统计和置换检验串成了一条完整流程适合手里已有标准空间脑图、想验证基因贡献的研究生也适合工程岗评估自己的数据能不能和转录组数据对上。核心不难坑全在坐标和批次效应上下面逐层拆开。2. 艾伦人类大脑图谱的数据结构六个供体、五万探针与MNI152坐标2.1 微阵列表达矩阵与样本注释先分清行、列和坐标列AHBA微阵列数据来自六个成年捐献者大脑每个大脑的左、右半球各取几百个样本点位置覆盖皮层和皮层下结构所以六个人加起来的样本量通常在两千到三千出头。表达量用Affymetrix Human Genome U133 Plus 2.0芯片测量探针数五万多个。下载包里最核心的三个文件是表达矩阵、探针注释、样本注释。% 三个核心文件一次性读入 probeInfo readtable(Probes.csv); % 探针注释gene_symbol、探针ID等 sampleAnnot readtable(SampleAnnot.csv); % 样本注释donor、结构名、MNI坐标 expr readmatrix(MicroarrayExpression.csv); % 探针 × 样本 的表达矩阵 fprintf(表达矩阵: %d 探针 × %d 样本\n, size(expr,1), size(expr,2)); fprintf(样本注释列名: %s\n, strjoin(sampleAnnot.Properties.VariableNames, , ));这里有几件事必须确认。readmatrix对纯数值CSV最省事但AHBA官网给的表达矩阵第一列偶尔带着探针ID文本直接读会得到NaN行。我一般先readtable看一眼前几行再决定或者显式指定readmatrix(..., NumHeaderLines, 1)。更关键的是SampleAnnot.csv里的空间信息坐标列名通常是mri_x、mri_y、mri_z单位是毫米。有人会把坐标和体素索引混为一谈后面所有采样全错这个问题在第4章单独说。2.2 从原始CSV到MATLAB结构体读取、清洗和探针过滤原始表达矩阵里有些样本的RNA质量差芯片检出率低直接拿进去算相关会引入噪声。AHBA同时提供了一个叫PACall.csv的文件记录每个探针在每个样本里是否显著高于背景。我一般会用这个做一层硬过滤检测不到的表达直接置为NaN不参与后续相关计算。% PACall 是present/absent call矩阵与expr同维度 pac readmatrix(PACall.csv); expr(pac 0.5) NaN; % 低于检测阈值的样本点置为缺失 % 统一放进一个结构体方便后续函数直接传 ahba struct(); ahba.expr single(expr); % 单精度内存直接减半 ahba.probeInfo probeInfo; ahba.sampleAnnot sampleAnnot; ahba.donorID sampleAnnot.donor; ahba.coords_mm [sampleAnnot.mri_x, sampleAnnot.mri_y, sampleAnnot.mri_z]; ahba.structure sampleAnnot.structure_name; save(ahba_clean.mat, ahba, -v7.3); % 超过2GB用v7.3格式single()这一步在五万探针级别的数据上非常值得做double 矩阵占的内存正好翻倍很多笔记本就是在这爆的内存。-v7.3的意义在于突破传统MAT格式的4GB上限数据量大时保存失败多半是因为忘了这个开关。清洗时还有一个常见误用有人按基因名直接删掉表达全为NaN的探针这没问题但千万别删样本注释行——样本注释和表达矩阵的列一一对应删错一行坐标就全错位了。2.3 坐标体系是前提MNI152毫米坐标怎么换算成体素索引AHBA样本坐标在MNI152空间而你下载的MRI图大概率也是NIfTI格式但NIfTI文件保存的是体素网格索引从1开始体素尺寸有1mm、2mm、3mm之分方向余弦也可能不一样。把毫米坐标直接除以体素大小当成索引是新手最容易犯的错。正确的做法是用NIfTI头里的齐次变换矩阵做世界坐标到体素坐标的逆变换。function [ix, iy, iz] mni2voxel(coords_mm, niiInfo) % 用nifti头里的Transform把MNI毫米坐标转成体素索引 % niiInfo.Transform是体素-世界的4x4齐次矩阵 T niiInfo.Transform.T; n size(coords_mm, 1); world [coords_mm, ones(n,1)]; vox T \ world; % 世界-体素左除等价于矩阵求逆 ix round(vox(1,:)); iy round(vox(2,:)); iz round(vox(3,:)); % 越界索引置NaN调用处直接跳过 out ix 1 | ix niiInfo.ImageSize(1) | ... iy 1 | iy niiInfo.ImageSize(2) | ... iz 1 | iz niiInfo.ImageSize(3); ix(out) NaN; iy(out) NaN; iz(out) NaN; endniiInfo.Transform.T是MATLAB从NIfTI头里解析出来的体素到世界变换矩阵行列顺序和SPM里spm_matrix的约定略有差异但本质一样。要用T \ world而不是T * world因为方向反了T 是把体素索引映射到毫米坐标现在要反着来。这个细节我见过不止一次被写反结果所有脑区坐标偏移了好几个厘米。3. 关联分析主流程球面采样、探针聚合、空间置换检验3.1 把MRI体素值插值到AHBA样本坐标上球面平均采样AHBA芯片检测的组织块不是单个体素而是毫米级别的立体取样核团。对应地从MRI图取值也应该以样本坐标为中心取一个小半径球体内的体素平均而不是只取最近的那个体素。尤其是2mm或3mm的MRI图单点取值会带来很大的采样误差。img niftiread(cortical_thickness.nii); info niftiinfo(cortical_thickness.nii); [ix, iy, iz] mni2voxel(ahba.coords_mm, info); % 球面半径取1.5mm换算成体素数 radiusVox max(1, round(1.5 / info.PixelDimensions(1))); nSamples numel(ix); mriVals nan(nSamples, 1); for s 1:nSamples if isnan(ix(s)), continue; end xr max(1, ix(s)-radiusVox) : min(info.ImageSize(1), ix(s)radiusVox); yr max(1, iy(s)-radiusVox) : min(info.ImageSize(2), iy(s)radiusVox); zr max(1, iz(s)-radiusVox) : min(info.ImageSize(3), iz(s)radiusVox); patch img(xr, yr, zr); % 去掉背景0值再取平均 vals patch(patch 0); if ~isempty(vals) mriVals(s) mean(vals(:)); end end ahba.mriVals mriVals;这里一个隐藏坑背景体素在很多MRI图里是0但脑脊液在某些序列里可能也是低值甚至接近0。严格做法是再套一个灰质掩膜只平均落在灰质内的体素否则白质和脑室的低信号会把皮层样本值往下拽。如果拿到的图是任务激活z-score图0值本身有意义不能随便当背景剔除——先确认图的类型再决定过滤策略。3.2 探针到基因的聚合策略选方差最大还是平均值芯片上一个基因往往对应多个探针组同一个基因的表达谱在不同探针下可能不完全一致。要做基因层面的关联分析必须先决定怎么把多探针合并成一个值。常见做法有两种取平均或选跨样本方差最大的探针。前者稳健后者信噪比高但选最大方差探针时有数据窥探的风险。我一般优先用最大方差探针因为它在空间上最能体现该基因的真实差异。geneSymbols probeInfo.gene_symbol; [geneList, ~, geneIdx] unique(geneSymbols); nGenes numel(geneList); exprGene zeros(nGenes, size(ahba.expr, 2), single); for g 1:nGenes rows find(geneIdx g); if numel(rows) 1 exprGene(g, :) ahba.expr(rows, :); else % 标准差最大的探针代表该基因 sd std(ahba.expr(rows, :), 0, 2, omitnan); [~, imax] max(sd); exprGene(g, :) ahba.expr(rows(imax), :); end end ahba.exprGene exprGene; ahba.geneNames geneList;这段逻辑里有个细节std的omitnan选项在R2019a之后才稳定可用。如果探针在部分样本里是NaN而这一探针恰好方差最大选出来的代表探针可能本身缺了一大半数据导致后续相关分析的有效样本数骤降。遇到这种情况我会退一步改取非NaN样本数最多的那个探针再比较方差。总之聚合策略要写进方法里审稿人一定会问。3.3 皮尔逊相关的统计推断普通p值、FDR和空间置换检验到了这一步数据是干净的MRI值、基因表达值、样本坐标都对齐了。核心计算其实就是逐基因算皮尔逊相关。但这里有两个统计陷阱必须处理多重比较校正以及空间自相关导致的自由度虚高。nGenes size(exprGene, 1); rAll nan(nGenes, 1); pAll nan(nGenes, 1); valid ~isnan(mriVals) sum(~isnan(exprGene), 1) 30; for g 1:nGenes y exprGene(g, :); ok valid ~isnan(y); if sum(ok) 30, continue; end [r, p] corr(mriVals(ok), y(ok)); rAll(g) r; pAll(g) p; end % FDR校正用生物信息学工具箱的mafdr或者下载fdr_bh adjP mafdr(pAll, BHFDR, true);普通p值计算非常容易但解释时要格外小心。AHBA样本在空间上高度相关相邻样本点的MRI值和基因表达都不独立所以经典t检验给出的自由度是虚高的p值会偏小。更稳的做法是拿置换检验兜底随机打乱样本坐标重算相关得到零分布再求经验p值。rObs rAll(targetGeneIdx); nPerm 5000; rNull zeros(nPerm, 1); for k 1:nPerm permIdx randperm(nSamples); rNull(k) corr(mriVals(valid), exprGene(targetGeneIdx, valid)); end pSpin (sum(abs(rNull) abs(rObs)) 1) / (nPerm 1); fprintf(观察r%.3f, 空间置换p%.4f\n, rObs, pSpin);置换检验这里有个边界条件打乱的是样本标签还是空间坐标结论不同。常见做法是打乱坐标与表达值的配对关系保留样本坐标间的空间相关结构这样零分布能反映空间自相关的影响。5000次置换在2万基因上跑不现实所以只对FDR筛选后的少数目标基因做复核别对所有基因跑。3.4 逐供体分析与元分析处理批次效应的标准做法六个供体的数据合在一起算相关问题在于供体间的全局表达差异死亡到取材的时间、RNA质量、芯片批次都会让表达谱整体偏移。这种批次效应很容易造成假阳性。标准做法是先按供体分别算相关再用Fisher z变换做固定效应元分析组合六个独立估计。donors unique(ahba.donorID); nd numel(donors); zVec zeros(nd, 1); wVec zeros(nd, 1); for d 1:nd sel strcmp(ahba.donorID, donors{d}); ok valid sel; if sum(ok) 30, continue; end [rd, ~] corr(mriVals(ok), exprGene(targetGeneIdx, ok)); zVec(d) atanh(rd); % Fisher z变换 wVec(d) sum(ok) - 3; % 每个供体的近似权重 end zMeta sum(zVec .* wVec) / sqrt(sum(wVec.^2)); pMeta 2 * (1 - normcdf(abs(zMeta))); fprintf(逐供体元分析: z%.3f, p%.3e\n, zMeta, pMeta);权重wVec用的是“有效样本数减3”这种近似严格讲应该用组内方差的倒数但对微阵列表达数据来说这个近似已经能控制住供体层面的异质性。逐供体分析还有个额外价值能看出某个显著结果是不是被单个供体主导。如果六个供体里只有一个人贡献了全部显著性这个结果基本不可复现趁早放弃。4. 避坑手册坐标错位、探针冗余、批次效应的五个真实排错记录4.1 坐标错位所有提取值都是背景现象做完球面采样后mriVals里超过六成是NaN连额叶样本都取不到有效值。原因MRI图文件名里带MNI152但实际上是从别的模板配准来的方向余弦和原点对不上另一种更常见的写法是把MNI毫米坐标直接除以体素大小忽略了NIfTI头的平移量。解决先用niftiinfo打印Transform手工验证几个已知脑区的坐标换算是否和FSLeyes一致如果明显错位就把MRI图先配准到ICBM152非对称2mm模板再用配准后的图重新采样。配准用SPM12的spm_normalise或者FSL的flirt都行关键是配准目标模板要和AHBA坐标所在的MNI152空间一致。4.2 同一基因多个探针相关结果互相矛盾现象同一个基因有五个探针五个r值从-0.3到正0.4都有结论悬空。原因不同探针覆盖同一转录本的不同外显子区域部分探针存在交叉杂交表达谱并不一致盲目取平均会把真实信号稀释掉。解决默认采用跨样本标准差最大的探针并在方法部分写明如果最大方差探针的有效样本数太少就改用非NaN样本最多的探针。不要看哪个探针p值小就留哪个那是典型的数据窥探审稿人一眼就能看出来。4.3 合并六个供体后相关一片红可能是批次效应在作怪现象把所有样本合并计算结果几乎每个基因都显著r值普遍在0.5以上脑图上的模式千篇一律。原因供体间的全局表达差差异直接混进了相关里。六个供体在表达均值上的差异如果和MRI值的空间分布有偶然重叠就会产生大片伪显著。解决先把表达数据在每个供体内部做z-score标准化再合并计算更稳妥的是直接走 3.4 节的逐供体元分析流程让每个供体的贡献单独估计。结果报告里至少标注“已按供体标准化”或“逐供体元分析”否则复现的人第一反应就是批次效应没处理。4.4 置换检验和普通p值结论相反现象普通相关p值小于0.001FDR也通过了但空间置换检验p值到了0.07逐供体元分析更是不显著。原因样本空间自相关让有效自由度远小于样本数经典t检验的p值系统性偏小置换检验保留了空间结构给出的零分布更贴近实际。解决以空间置换检验或逐供体元分析的结果为准普通p值只当作筛选候选基因的粗指标。写文章时把两种方法的结果都放出来审稿人反而更信任。数值上置换次数低于2000时p值分辨率不够至少跑到5000。4.5 内存不足、矩阵太大现象readmatrix加载表达矩阵时MATLAB直接转圈十几分钟或者save时报内存不足。原因默认double精度下五万探针乘两三千样本就有上GB占用量某些老版本MAT文件格式有4GB限制保存大矩阵会直接失败。解决读取后立刻转single精度表达矩阵按供体拆分处理不要总想着一次全load进内存保存用save(..., -v7.3)。如果只是验证少数基因读取时用readmatrix的Range参数只加载目标行能省掉大量等待时间。5. 进阶技巧用标签映射快速生成基因-脑区共定位报告关联分析做完通常要回答“这个基因主要在哪些脑区表达”。与其逐个样本点看坐标不如把AHBA样本映射到解剖图谱标签上按脑区聚合出表达模式。这里推荐一个简单可靠的做法把AHBA样本坐标投影到Desikan-Killiany图谱或AAL3图谱的NIfTI文件上用图谱标签给每个样本点打上脑区名再按脑区统计目标基因的平均表达。atlas niftiread(DK90_2mm.nii); atlasInfo niftiinfo(DK90_2mm.nii); [ax, ay, az] mni2voxel(ahba.coords_mm, atlasInfo); nSamples numel(ax); roiLabel zeros(nSamples, 1); for s 1:nSamples if isnan(ax(s)), continue; end roiLabel(s) atlas(ax(s), ay(s), az(s)); end % 按脑区聚合目标基因表达 geneIdx find(strcmp(ahba.geneNames, COX7A2)); valid ~isnan(roiLabel) ~isnan(exprGene(geneIdx, :)); roiData table(ahba.structure(valid), roiLabel(valid), exprGene(geneIdx, valid)); roiData.Properties.VariableNames {structure, label, expression}; roiSummary groupsummary(roiData, structure, mean, expression);这段代码用了groupsummaryMATLAB R2018a之后都有。如果不想依赖图谱配准误差也可以直接用ahba.structure里的解剖名聚合但结构名层级多、命名不统一用图谱标签反而更干净。这一套做完你手头就有了从MRI脑图到基因表达、再到脑区富集的三层结果。从那以后我每次跑关联之前都强制自己核对三个数样本坐标的质控、有效样本数、供体数。摘要里如果缺这三个数我基本不会签字。希望帮到你。本文还有配套的精品资源点击获取