本帖最后由 Ogiwara 于 2024-5-10 21:02 编辑
小白第一次发帖，请各位大佬多多海涵 去年11月好不容易注册一个号，不要我长时间不发帖又给我收回了
事情是这样的：lz上课拿到一个数据集，是很多NMR谱按照0.04ppm的区间进行切块做成的excel表，课题要求做PCA、PLS-DA分析、STOCSY分析并且分析出哪些区块是互相关联的，又对应的是哪些化合物。PCA、PLS-DA好使，网上都有现成的R语言工具，STOCSY也用了R语言把数据拿到手了。放进excel表一看，我去，怎么这么多格子，这一看就不是人能干的活啊！
STOCSY数据，好多格子
图1 逆天STOCSY表
于是lz我就打开了python开始写脚本了。
既然要分析出哪些区块是互相关联的，那么就要把那些值不是0的格子全部读取出来，再读取对应的首行和首列格子里面的值，之后分别存到表格里面，再把表格里面的值分别填到另外一份工作表中。lz之前写过群机器人的lua脚本，所以还是看得懂代码的，就是写出来的东西不一定能跑。python的xlwings模块也没用过，于是就打开百度上网抄了几行代码：
[Python] 纯文本查看 复制代码 import xlwings as xw
app = xw.App(visible=True,add_book=False)
#不显示Excel消息框
app.display_alerts=False
#关闭屏幕更新,可加快宏的执行速度
app.screen_updating=False
wb = app.books.open("D:\\Study\\NMR\\Prac_3\\correlation.xlsx")
# 输出打开的excle的绝对路径
cort = wb.sheets["correlation_table"]
opposi = wb.sheets["outputposi"]
opnega = wb.sheets["outputnega"]
#读取三个表，其中cort是数据源，posi放在0~1区间的一对化学位移，nega放小于0的一对化学位移。
#接下来首先遍历表中B3:IC238来寻找符合条件的值，寻找到之后就把对应值那一行第一个值放进arrow表内，那一列第一个值放进column表内。
opposi_row = []
opposi_col = []
opnega_row = []
opnega_col = []
pval_posi = []
pval_nega = []
#先建空的列表
print("loading data...")
#接下来是遍历模块
rows, cols = cort.used_range.shape
#先把有用的总行数、总列数全部读取出来
#B3:IC238应该是从第3行，第2列开始读取
for row in range(3, rows):
for col in range(2, cols):
pval = cort.range(row,col).value
if pval < 0:
cs_row = cort.range(row,1).value
cs_col = cort.range(1,col).value
#表中值为负数时就把这个先赋值一个变量
opnega_row.append(cs_row)
opnega_col.append(cs_col)
pval_nega.append(pval)
#再把变量塞进表里
else:
if pval > 0 and pval < 1:
#正数也是一样
cs_row = cort.range(row,1).value
cs_col = cort.range(1,col).value
opposi_row.append(cs_row)
opposi_col.append(cs_col)
pval_posi.append(pval)
opposi.range('A1').value = opposi_row
opposi.range('A2').value = opposi_col
opposi.range('A3').value = pval_posi
opnega.range('A1').value = opnega_row
opnega.range('A2').value = opnega_col
opnega.range('A3').value = pval_nega
#把列表里面的东西全部放进输出里面
print("done")
wb.save()
wb.close()
# 退出excel程序
app.quit()
第一次用xlwings居然没有报错，跑了十几分钟打开表格一看都给我处理完了。虽然都是一行行存的，但是之后复制转置一下就好了。
拿到的数据
图2 脚本输出数据，第一二列是化学位移，第三列是对应相关系数
但是再看了一下数据几千行还是有点多，而且没有峰的地方仪器噪声硬是给我关联上了，再说了课题要求是把化合物之间的联系全部列出来，所以进一步的数据处理还是必要的。于是lz又打开老教授给的NMR谱文件开始了艰苦且感觉并不准确的手动化合物分配，拿到了一张化学位移表。
据说R语言里面有一个叫Batman的组件可以自动分配NMR谱里面的化合物，可是lz并没有安装R语言的环境，只能用试用版软件手动分配，花了很多没用的时间，再也不想做第二次了
之后又写了个脚本，具体思路是读图2输出的数据的第一列，并在同时在化学位移表里面遍历找为第一列数据+-0.04ppm的峰，找到就读取对应第一列的化合物名字存变量，同时读取第二列数据，并且再在这个基础上再在化学位移表里面遍历拿化合物名称存另一个变量，如果两个化合物名称都读取到了就把p值也存进变量里，然后统统塞表格再放到新工作表里面。以上面抄到的代码为蓝本，lz又修修补补自己的代码，写出了惊人的5+4个for嵌套循环：
[Python] 纯文本查看 复制代码 import xlwings as xw
app = xw.App(visible=True,add_book=False)
#不显示Excel消息框
app.display_alerts=False
#关闭屏幕更新,可加快宏的执行速度
app.screen_updating=False
wb = app.books.open("D:\\Study\\NMR\\Prac_3\\correlation.xlsx")
# 输出打开的excle的绝对路径
opposi = wb.sheets["outputposi"]
opnega = wb.sheets["outputnega"]
metaref = wb.sheets["meta_shift"]
metaop = wb.sheets["meta_cor"]
#读取四个表，opposi和opnega都是上次输出的有联系的区块，但是都是0.04ppm这么切的，之后也只要0.04设为阈值就行。
#接下来遍历metaref，opposi和opnega，先遍历opposi和opnega的第一列，找到一样的之后把对应化合物塞入metacor1，之后定义值为第二列的化学位移然后到metaref里面找有没有在阈值内的，找到就把对应化合物塞入metacor2，p值塞入metapval。
metacor1 = []
metacor2 = []
metapval = []
#先把空表格弄好
print("loading data...")
#然后开始迭代
rows, cols = metaref.used_range.shape
posirows, posicols = opposi.used_range.shape
negarows, negacols = opnega.used_range.shape
#先把有用的总行数、总列数全部读取出来
for row in range(1, rows):
for col in range(2, cols):
metacs = metaref.range(row,col).value
#读取化学位移表中代谢物化学位移
for posirow in range (1,posirows):
posics = opposi.range(posirow,1).value
#读取正表中第一列化学位移
compareval = metacs - posics
if compareval <= 0.04 and compareval >= -0.04:
#如果条件成立
metaname = metaref.range(row,1).value
#把化合物名字存在变量里
paircs = opposi.range(posirow,2).value
#配对的区域的化学位移值也存在变量里
metap = opposi.range(posirow,3).value
#对应的p值也存一下
for rowp in range(1,rows):
for colp in range(2,cols):
#回到对照表继续遍历
metapaircs = metaref.range(rowp,colp).value
#再读取一遍化学位移值
comval = paircs - metapaircs
if comval <= 0.04 and comval >= -0.04:
#如果符合条件
metapairname = metaref.range(rowp,1).value
#把配对的东西存变量
metacor1.append(metaname)
metacor2.append(metapairname)
metapval.append(metap)
for negarow in range (1, negarows):
negacs = opnega.range(negarow,1).value
comp = metacs - negacs
if comp <= 0.04 and comp >= -0.04:
metaname = metaref.range(row,1).value
#把化合物名字存在变量里
paircs = opnega.range(negarow,2).value
#配对的区域的化学位移值也存在变量里
metap = opnega.range(negarow,3).value
#对应的p值也存一下
for rowp in range(1,rows):
for colp in range(2,cols):
#回到对照表继续遍历
metapaircs = metaref.range(rowp,colp).value
#再读取一遍化学位移值
comval = paircs - metapaircs
if comval <= 0.04 and comval >= -0.04:
#如果符合条件
metapairname = metaref.range(rowp,1).value
#把配对的东西存变量
metacor1.append(metaname)
metacor2.append(metapairname)
metapval.append(metap)
metaop.range('A1').value = metacor1
metaop.range('A2').value = metacor2
metaop.range('A3').value = metapval
#变量存一下
print("done")
wb.save()
wb.close()
# 退出excel程序
app.quit()
这个代码加上几千行的数据，导致lz的代码虽然没有出bug但是跑这个代码的破笔记本被硬控了两个小时，而代码写完已经是半夜0点，于是lz开着笔记本先睡了一觉，醒来电脑已经帮我处理好了：
化合物之间关联
图3 化合物配对与对应系数
这样看着应该是好了，但数据还是太多，而且排除了值为1的峰（就是峰和峰本身的变化系数）还有重复的，有很多是因为化合物本身有自旋耦合出了裂分峰，还有是因为化合物有好几簇峰对应。而且这些东西里面一些化合物之间的联系还是有正有负的，于是lz起了把同名化合物对应的变化系数全部求平均然后塞二维表格的想法。
于是又编了一个脚本，想法是搞个嵌套字典，第一列的化合物名称塞进外面字典的key值里面，第二列的化合物塞进里面的字典的key里面，然后建个空表一点点塞东西。之后给字典里面的列表求平均，然后遍历放结果的表，列符合外面字典的key，行符合里面字典的key就把对应的平均值塞进去。花了半个小时编了之后发现表格里面没有东西，又花了一个小时才发现是因为后面验证条件的工作表写错了，再花了一个小时解决了每次读新化合物都会把里面字典清空的问题，最终可以正常运行的版本如下：
[Python] 纯文本查看 复制代码 import xlwings as xw
import numpy as np
app = xw.App(visible=True,add_book=False)
#不显示Excel消息框
app.display_alerts=False
#关闭屏幕更新,可加快宏的执行速度
app.screen_updating=False
wb = app.books.open("D:\\Study\\NMR\\Prac_3\\correlation.xlsx")
# 输出打开的excle的绝对路径
metacor = wb.sheets["meta_cor"]
metares = wb.sheets["cor_result"]
#读取工作表
#读取metacor的第一列同时读取metacor的第二列、第三列。搞一个字典，按照第一、第二列名称来存储第三列的p值。之后遍历字典按照第一列寻找行，第二列寻找列，然后把对应的表内的p值平均值塞进去。
database = {}
#建一个字典
rows, cols = metacor.used_range.shape
#读取一下metacor的行列数
print("reading metacor...")
for row in range(2, rows):
metaname1 = metacor.range(row,1).value
metaname2 = metacor.range(row,2).value
metapval = metacor.range(row,3).value
if metaname1 not in database:
database[metaname1] = {metaname2: []}
database[metaname1][metaname2].append(metapval)
else:
dict_ = database[metaname1]
if metaname2 not in dict_:
database[metaname1][metaname2] = []
database[metaname1][metaname2].append(metapval)
else:
database[metaname1][metaname2].append(metapval)
#把所有值塞进字典
resrows, rescols = metares.used_range.shape
print("inserting metares...")
for key in database.keys():
dict_ = database[key]
for key1 in dict_.keys():
pvalave = np.mean(database[key][key1])
for resrow in range(2, resrows):
compound1 = metares.range(resrow,1).value
if f"{key}" == compound1:
for rescol in range (2, rescols):
compound2 = metares.range(1,rescol).value
if f"{key1}" == compound2:
metares.range(resrow, rescol).value = pvalave
print("done")
wb.save()
wb.close()
# 退出excel程序
app.quit()
平常地跑了出来。
设置了一下条件格式，但是感觉数据还是有点多。然后又写了个脚本把工作表里面绝对值小于0.7的值全部改成0，还把自己和自己对应的系数也全改成0了，这下清楚多了。
以下是脚本：
[Python] 纯文本查看 复制代码 import xlwings as xw
app = xw.App(visible=True,add_book=False)
#不显示Excel消息框
app.display_alerts=False
#关闭屏幕更新,可加快宏的执行速度
app.screen_updating=False
wb = app.books.open("D:\\Study\\NMR\\Prac_3\\correlation.xlsx")
# 输出打开的excle的绝对路径
metapro = wb.sheets["cor_result_process"]
rows, cols = metapro.used_range.shape
print("processing...")
for row in range(2,rows):
for col in range(2,cols):
pval = metapro.range(row,col).value
if pval < 0.7 and pval > -0.7:
metapro.range(row,col).value = 0
print("done")
wb.save()
wb.close()
# 退出excel程序
app.quit()
总之就把数据处理完了。一些地方的处理可能存在问题，比如我分配化合物峰的时候直接填的是中心的化学位移，而没有考虑到裂分峰中心可能不会有那么大的变化，不过管他呢，我又不发文章，只是完成作业而已 继续写论文去了 | 
发表于 2024-5-10 19:35
发表于 2024-5-10 20:38
