前面我们给大家分享了一个超详细的python文献代码：（IF=50.0/Q1）顶刊杂志的python代码复现资源分享 ，今天就来学习它！势必要拿下！想进python群一起学习交流的加我微信：Biotree123。文献标题为：《Regenerative lineages and immune-mediated pruning in lung cancer metastasis》。
回顾 单细胞数据上传到了GEO数据库：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE123902
注释后的文件：
https://s3.amazonaws.com/dp-lab-data-public/lung-development-cancer-progression/PATIENT_LUNG_ADENOCARCINOMA_ANNOTATED.h5
https://s3.amazonaws.com/dp-lab-data-public/lung-development-cancer-progression/MOUSE_LUNG_ADENOCARCINOMA_METASTASIS_ANNOTATED.h5
代码下载的地方：https://github.com/dpeerlab/lung-development-cancer-progression
如有文件和代码下载困难，可加微信：Biotree123来获取~
第一篇：（IF=50.0/Q1）顶刊python代码学习（一）：数据背景，环境配置与数据加载
Note：
下面一些代码会比较长，代码里面给了一个 ipynb 的文件，直接在上面运行即可，不需要复制这里面的。 我们这里简单做一些代码解读，还有运行不成功的代码更新，确保都能正常运行。
今天的学习内容： 附图2中的d/e/h. h图主要对细胞进行初始大类注释，所有细胞最初被分配到三个大细胞类别之一：淋巴细胞、髓样细胞和上皮细胞/基质细胞类型。这种分配是根据细胞对经典免疫细胞、髓样细胞、间充质细胞和上皮细胞标记物的表达情况进行的。
1）加载数据 见稿子：（IF=50.0/Q1）顶刊python代码学习（一）：数据背景，环境配置与数据加载
2）评估库参数 设置绘图参数：
# exec 是 Python 的内置函数，用于执行存储在字符串或对象中的 Python 代码。 exec( 'QUERY = DF_{}' .format(subset_type)) plt.figure(figsize = ( 9 , 3 )) # 设置图形的宽度为 9 英寸，高度为 3 gs1 = gridspec.GridSpec( 1 , 3 ) # 创建一个 1 行 3 列的网格布局，即总共可以放置 3 个子图。 gs1.update(wspace= 0.5
, hspace= 0.7 ) # set the spacing between axes. 设置子图之间的水平间距为 0.5，垂直间距为 0.7。 子图1：注意这里的代码更新 在 Python 的 pandas 库中， DataFrame 对象并没有 select 方法。你可能想使用的是 loc 或 iloc 方法来选择数据。 select 方法是 pandas 早期版本中的方法，但在较新的版本中已经被废弃。
# (1) PLOT LOG10 LIBRARY SIZE ax = plt.subplot(gs1[0]) bins = np.linspace(np.log10(QUERY.values.sum(axis=1).min()), np.log10(QUERY.values.sum(axis=1).max())*0.95, 100) for ind, label in enumerate(np.unique(QUERY.index.get_level_values( 'Legend' ))): # 选择对应于 'Legend' 级别为 label 的行 # x = np.log10(QUERY.select(lambda x: x[1] in label).values.sum(axis=1)) # log cell size, 原来的代码，修改为下面的 x = np.log10(QUERY.loc[QUERY.index.get_level_values( 'Legend' ) == label].values.sum(axis=1)) # log cell size ，修改后的 # plt.hist(x, bins, alpha=0.5, label=label,color = np.divide(tuple(hex(FLATUI_SAMPLE[ind]).rgb),255)) # 这里缺少一个colors的库，网络上找不到，就不设置颜色了 plt.hist(x, bins, alpha=0.5, label=label) plt.xticks(rotation=70) plt.ylabel( 'Frequency' ) plt.xlabel( '\nLog10 Cell Size' ) plt.grid(True) sns.despine() 子图2： # (2) PLOT LOG NUM UNIQUE GENES ax = plt.subplot(gs1[ 1 ]) bins = np.linspace(np.log10((np.sum(QUERY > 0 ,axis= 1 )).min()), np.log10((np.sum(QUERY > 0 ,axis= 1 )).max())* 0.95 , 100 ) for ind, label in enumerate(np.unique(QUERY.index.get_level_values( 'Legend' ))): # x = np.log10(np.sum(QUERY.select(lambda x: x[1] in label) > 0,axis=1)) # plt.hist(x, bins, alpha=0.5, label=label,color = np.divide(tuple(hex(FLATUI_SAMPLE[ind]).rgb),255)) subset = QUERY.loc[QUERY.index.get_level_values( 'Legend' ) == label] x = np.log10(np.sum(subset > 0 , axis= 1 )) # log number of unique genes plt.hist(x, bins, alpha= 0.5 , label=label) plt.xticks(rotation= 70 ) plt.ylabel( 'Frequency' ) plt.xlabel( '\nLog10 Unique genes' ) plt.grid( True ) sns.despine() 子图3： # (3) PLOT LOG NUMBER OF CELLS CONTRIBUTING TO EACH GENE ax = plt.subplot(gs1[2]) data = np.log10(np.sum(QUERY.values > 0,axis=0)) bins = np.linspace(0, data.max()*0.95, 20) for ind, label in enumerate(np.unique(QUERY.index.get_level_values( 'Legend' ))): # x = np.log10(np.sum(QUERY.select(lambda x: x[1] in label) > 0,axis=0)) # x[(np.isinf(x)) | (np.isnan(x))] = 0 # plt.hist(x, bins, alpha=0.5, label=label,color = np.divide(tuple(hex(FLATUI_SAMPLE[ind]).rgb),255)) # 选择对应于 'Legend' 级别为 label 的行 subset = QUERY.loc[QUERY.index.get_level_values( 'Legend' ) == label] x = np.log10(np.sum(subset.values > 0, axis=0)) x[(np.isinf(x)) | (np.isnan(x))] = 0 # 处理无穷大和NaN值 plt.hist(x, bins, alpha=0.5, label=label) plt.xticks(rotation=70) plt.ylabel( 'Frequency' ) plt.xlabel( '\nLog. Cells/Gene' ) plt.grid(True) sns.despine() 三张图的代码放在一起运行得到的结果拼图：
这三个图分别展示了每个细胞中总umi数，每个细胞中表达的基因数，前两者的比值。不同颜色表示不同样本。
3）初始大类注释 1.读取细胞marker 数据在 CELL_TYPES_WANSLEEBEN_HOGAN_2013.csv 表格，作者在提供的代码data文件夹里面有：
# 从一个 Excel 文件中加载标准的细胞类型标记基因。 # 这个文件对应于 SUPPLEMENTARY_TABLE2/TAB1 (Gene Signatures from ED Fig. 1i) genesets = pd.read_csv( './data/CELL_TYPES_WANSLEEBEN_HOGAN_2013.csv' ,header= 'infer' ) # 让 pandas 自动推断文件的头部行（通常是第一行）作为列名 print(genesets.shape[ 1 ])
# 输出多少列 genesets 2.生成每个细胞类型基因集中所有基因表达总和的打分矩阵： # ROW INDEX meta = 'META_CELL_TYPE' FLATUI_PLOT = FLATUI_METACELL title = os.path.basename(path_to_genesets)[: -4 ] # 从 genesets 数据框中提取所有列名，这些列名对应不同的细胞类型 title # QUERY 为INDF_ALL # GENERATE MATRIX OF CUMULATIVE EXPRESION PER CELL TYPE exec( 'QUERY = INDF_{}' .format(subset_type)) titles = list(genesets.columns) #['EPITHELIAL','MESENCHYAL','CD45'] tmp = pd.DataFrame(index = QUERY.index) # 初始化一个空的 DataFrame，其索引与 QUERY 相同 # 循环：对于每个细胞类型（title），从 genesets 中提取对应的基因列表，过滤掉 NaN 值，并找出这些基因在 QUERY 中的交集。然后计算这些基因在每个细胞中的累积表达值，并将结果存储在 tmp 数据框中。 for title in titles: genes = genesets.values[:,np.where(genesets.columns==title)[ 0 ][ 0 ]] # genes x pathways genes = [x for x in genes if str(x) != 'nan' ] detected_genes = list(set(genes).intersection(set(QUERY.columns))) vals = QUERY[detected_genes] exec( 'tmp[\'{}_SCORE\'] = (np.nansum(vals,axis=1))' .format(title)) genes = list(tmp.columns) 第三步的颜色设置去掉了。
4. 构建热图数据 对上面的基因集累积表达值进行了归一化处理：进行标准化（Z-score）
# CONSTRUCT HEATMAP DATA heatmap_data = pd.DataFrame(data = zscore(tmp.values,axis=0),columns = genes) # yticks 和 xticks：分别存储热图的行和列标签。 yticks = heatmap_data.index xticks = heatmap_data.columns heatmap_data 5.对上面的归一化值进行行列聚类 # LINKAGE method = 'average' # average, single metric = 'cosine' # cosine linkage = hc.linkage(heatmap_data, method=method, metric = metric) row_linkage = deepcopy(linkage) linkage = hc.linkage(heatmap_data.T, method=method, metric = metric) col_linkage = deepcopy(linkage) # 对热图按照上面的聚类结果重新进行排序 r1 = hc.leaves_list(row_linkage) cl = hc.leaves_list(col_linkage) mat = heatmap_data.iloc[r1,cl] mat
行列已经发生了变化。
6.绘制聚类标记热图 颜色设置的地方注释掉了。
# VIEW CLUSTERED LABELED HEATMAP AND DENDROGRAM fig = plt.figure(figsize=( 18 , 10 )) # ADD ROW COLOR INDEX (CELL OF ORIGIN) ax1 = fig.add_axes([ 0.088 , 0.1 , 0.01 , 0.6 ]) # [x0,y0,width,height] x = 0 y = 0 # for c in row_colors[r1]: # pos = (x, y / len(r1)) # ax1.add_patch(patches.Rectangle(pos, 1, 1 / len(r1), color=c)) # if y >= len(r1)-1: # x += 1 # y = 0 # else: # y += 1 # plt.axis('off') # ADD MATRIX WITH GENE NAMES axmatrix = fig.add_axes([ 0.1 , 0.1 , 0.7 , 0.6 ]) im = axmatrix.matshow(mat, aspect= 'auto' , origin= 'lower' , cmap=plt.cm.RdBu_r,vmin= -3 ,vmax= 3 ) labels = list(mat.columns) axmatrix.xaxis.set_ticks_position( 'bottom' ) axmatrix.set_yticklabels([]); xtick = plt.xticks(range(len(labels)), labels, rotation = 90 , fontsize = 14 ) # ADD DENDROGRAM sch.set_link_color_palette([ '#808080 ' , '#808080 ' , '#808080 ' , '#808080 ' , '#808080 ' , '#808080 ' , '#808080 ' ]) ax2 = fig.add_axes([ 0.01 , 0.1 , 0.08 , 0.6 ]) # [x0,y0,width,height] Z1 = sch.dendrogram(row_linkage, orientation= 'left' ,above_threshold_color= '#808080 ' ) ax2.set_xticks([]) ax2.set_yticks([]) plt.axis( 'off' ) # ADD COLORBAR axcolor = fig.add_axes([ 0.81 , 0.1 , 0.01 , 0.6 ]) cbar = plt.colorbar(im, cax=axcolor) cbar.ax.get_yaxis().set_ticks([])
figure_label = 'heatmaps/user_defined/_{}_{}_HEATMAP_{}_{}_{}' .format(meta,subset_type,title,method,metric) fn = FIG_output_stem + FN.replace( ".h5" , "" ) + figure_label d = os.path.dirname(fn) if not os.path.exists(d): os.makedirs(d) fig.savefig(fn + '.png' , bbox_inches= 'tight' ,dpi=fig_dpi) print(fn) 感觉没有必要这么复杂呀，经典的细胞类型很容易就区分出：淋巴细胞、髓样细胞和上皮细胞/基质细胞类型。
今天分享到这~
（未完待续）
