Draw 3D PCA and calculate background genes

gene_x 0 like s 1073 view s

Tags: plot, python

Download python3 script for generating the plot

Download raw data 1 for generating the plot

Download raw data 2 for generating the plot

PCA_3D_2.png

construct a data structure (merged_df) as above with data and pc

  1. library(ggplot2)
  2. data <- plotPCA(rld, intgroup=c("condition", "donor"), returnData=TRUE)
  3. write.csv(data, file="plotPCA_data.csv")
  4. #calculate all PCs including PC3 with the following codes
  5. library(genefilter)
  6. ntop <- 500
  7. rv <- rowVars(assay(rld))
  8. select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
  9. mat <- t( assay(rld)[select, ] )
  10. pc <- prcomp(mat)
  11. pc$x[,1:3]
  12. #df_pc <- data.frame(pc$x[,1:3])
  13. df_pc <- data.frame(pc$x)
  15. identical(rownames(data), rownames(df_pc)) #-->TRUE
  16. ## define the desired order of row names
  17. #desired_order <- rownames(data)
  18. ## sort the data frame by the desired order of row names
  19. #df <- df[match(desired_order, rownames(df_pc)), ]
  21. data$PC1 <- NULL
  22. data$PC2 <- NULL
  23. merged_df <- merge(data, df_pc, by = "row.names")
  24. #merged_df <- merged_df[, -1]
  25. row.names(merged_df) <- merged_df$Row.names
  26. merged_df$Row.names <- NULL  # remove the "name" column
  27. merged_df$name <- NULL
  28. merged_df <- merged_df[, c("PC1","PC2","PC3","PC4","PC5","PC6","PC7","PC8","PC9","PC10","PC11","PC12","PC13","PC14","PC15","PC16","PC17","PC18","PC19","PC20","PC21","PC22","PC23","PC24","PC25","PC26","PC27","PC28","group","condition","donor")]
  29. #results in 26PCs: merged_df <- merged_df[, c("PC1","PC2","PC3","PC4","PC5","PC6","PC7","PC8","PC9","PC10","PC11","PC12","PC13","PC14","PC15","PC16","PC17","PC18","PC19","PC20","PC21","PC22","PC23","PC24","PC25","PC26","group","condition","donor")]
  30. write.csv(merged_df, file="merged_df_28PCs.csv")
  31. #> summary(pc)  #--> variances of PC1, PC2 and PC3 are 0.6795 0.08596 0.06599
  32. #Importance of components:
  33. #                          PC1     PC2     PC3     PC4     PC5     PC6     PC7
  34. #Standard deviation     2.6011 0.92510 0.81059 0.67065 0.51952 0.46429 0.41632
  35. #Proportion of Variance 0.6795 0.08596 0.06599 0.04517 0.02711 0.02165 0.01741
  36. #Cumulative Proportion  0.6795 0.76548 0.83148 0.87665 0.90376 0.92541 0.94282
  37. #                           PC8     PC9    PC10    PC11    PC12    PC13    PC14
  38. #Standard deviation     0.38738 0.32048 0.27993 0.24977 0.20217 0.17316 0.16293
  39. #Proportion of Variance 0.01507 0.01032 0.00787 0.00627 0.00411 0.00301 0.00267
  40. #Cumulative Proportion  0.95789 0.96821 0.97608 0.98234 0.98645 0.98946 0.99213
  41. #                          PC15    PC16    PC17    PC18    PC19    PC20    PC21
  42. #Standard deviation     0.15058 0.13083 0.12449 0.08789 0.06933 0.06064 0.04999
  43. #Proportion of Variance 0.00228 0.00172 0.00156 0.00078 0.00048 0.00037 0.00025
  44. #Cumulative Proportion  0.99441 0.99612 0.99768 0.99846 0.99894 0.99931 0.99956
  45. #                          PC22    PC23    PC24    PC25     PC26
  46. #Standard deviation     0.04143 0.03876 0.03202 0.01054 0.005059
  47. #Proportion of Variance 0.00017 0.00015 0.00010 0.00001 0.000000
  48. #Cumulative Proportion  0.99973 0.99988 0.99999 1.00000 1.000000

draw 3D with merged_df using plot3D

  1. import plotly.graph_objects as go
  2. import pandas as pd
  3. from sklearn.decomposition import PCA
  4. import numpy as np
  6. # Read in data as a pandas dataframe
  7. #df = pd.DataFrame({
  8. #    'PC1': [-13.999925, -12.504291, -12.443057, -13.065235, -17.316215],
  9. #    'PC2': [-1.498823, -3.342411, -6.067055, -8.205809, 3.293993],
  10. #    'PC3': [-3.335085, 15.207755, -14.725450, 15.078469, -6.917358],
  11. #    'condition': ['GFP d3', 'GFP d3', 'GFP d8', 'GFP d8', 'GFP+mCh d9/12'],
  12. #    'donor': ['DI', 'DII', 'DI', 'DII', 'DI']
  13. #})
  14. df = pd.read_csv('merged_df_28PCs.csv', index_col=0, header=0)
  15. df['condition'] = df['condition'].replace("GFP+mCh d9/12", "ctrl LTtr+sT d9/12")
  16. df['condition'] = df['condition'].replace("sT+LTtr d9/12", "LTtr+sT d9/12")
  17. df['condition'] = df['condition'].replace("GFP d3", "ctrl LT/LTtr d3")
  18. df['condition'] = df['condition'].replace("GFP d8", "ctrl LT/LTtr d8")
  19. df['condition'] = df['condition'].replace("mCh d3", "ctrl sT d3")
  20. df['condition'] = df['condition'].replace("mCh d8", "ctrl sT d8")
  22. # Fit PCA model to reduce data dimensions to 3
  23. pca = PCA(n_components=3)
  24. pca.fit(df.iloc[:, :-3])
  25. X_reduced = pca.transform(df.iloc[:, :-3])
  27. # Add reduced data back to dataframe
  28. df['PC1'] = X_reduced[:, 0]
  29. df['PC2'] = X_reduced[:, 1]
  30. df['PC3'] = X_reduced[:, 2]
  32. # Create PCA plot with 3D scatter
  33. fig = go.Figure()
  35. #['circle', 'circle-open', 'square', 'square-open', 'diamond', 'diamond-open', 'cross', 'x']
  36. # if donor == 'DI' else marker=dict(size=2, opacity=0.8, color=condition_color, symbol=donor_symbol)
  38. #decrease diamond size to 6 while keep the circle as size 10 in the following code:
  39. #'rgb(128, 150, 128)'
  40. #I need three families of colors, always from light to deep, the first one should close to grey.
  41. #the first serie for 'ctrl LTtr+sT d9/12', 'LTtr+sT d9/12' 
  42. #the second serie for 'ctrl LT/LTtr d3', 'ctrl LT/LTtr d8', 'LT d3', 'LT d8', 'LTtr d3', 'LTtr d8'
  43. #the third serie for 'ctrl sT d3', 'ctrl sT d8', 'sT d3', 'sT d8', 'sT+LT d3'
  45. condition_color_map_untreated = {'untreated':'black'}
  46. donor_symbol_map_untreated = {'DI': 'circle-open', 'DII': 'diamond-open'}
  47. #condition_color_map = {'ctrl LTtr+sT d9/12': 'green', 'GFP d3': 'blue', 'GFP d8': 'red', 'GFP+mCh d9/12': 'green', 'LT d3': 'orange'}
  48. condition_color_map = {
  49.     'ctrl LTtr+sT d9/12': 'black',
  50.     'LTtr+sT d9/12': '#a14a1a',
  52.     'ctrl LT/LTtr d3': '#fdbf6f',
  53.     'ctrl LT/LTtr d8': '#ff7f00',
  54.     'LT d3': '#b2df8a',
  55.     'LT d8': '#33a02c',
  56.     'LTtr d3': '#a6cee3',
  57.     'LTtr d8': '#1f78b4',
  59.     'ctrl sT d3': 'rgb(200, 200, 200)',
  60.     'ctrl sT d8': 'rgb(100, 100, 100)',
  61.     'sT d3': '#fb9a99',
  62.     'sT d8': '#e31a1c',
  63.     'sT+LT d3': 'magenta'
  64. }
  65. donor_symbol_map = {'DI': 'circle', 'DII': 'diamond'}
  67. for donor, donor_symbol in donor_symbol_map_untreated.items():
  68.     for condition, condition_color in condition_color_map_untreated.items():
  69.         mask = (df['condition'] == condition) & (df['donor'] == donor)
  70.         fig.add_trace(go.Scatter3d(x=df.loc[mask, 'PC1'], y=df.loc[mask, 'PC2'], z=df.loc[mask, 'PC3'],
  71.                                    mode='markers',
  72.                                    name=f'{condition}' if donor == 'DI' else None,
  73.                                    legendgroup=f'{condition}',
  74.                                    showlegend=True if donor == 'DI' else False,
  75.                                    marker=dict(size=6 if donor_symbol in ['diamond-open'] else 10, opacity=0.8, color=condition_color, symbol=donor_symbol)))
  77. for donor, donor_symbol in donor_symbol_map.items():
  78.     for condition, condition_color in condition_color_map.items():
  79.         mask = (df['condition'] == condition) & (df['donor'] == donor)
  80.         fig.add_trace(go.Scatter3d(x=df.loc[mask, 'PC1'], y=df.loc[mask, 'PC2'], z=df.loc[mask, 'PC3'],
  81.                                    mode='markers',
  82.                                    name=f'{condition}' if donor == 'DI' else None,
  83.                                    legendgroup=f'{condition}',
  84.                                    showlegend=True if donor == 'DI' else False,
  85.                                    marker=dict(size=6 if donor_symbol in ['diamond'] else 10, opacity=0.8, color=condition_color, symbol=donor_symbol)))
  87. for donor, donor_symbol in donor_symbol_map.items():
  88.     fig.add_trace(go.Scatter3d(x=[None], y=[None], z=[None],
  89.                                mode='markers',
  90.                                name=donor,
  91.                                legendgroup=f'{donor}',
  92.                                showlegend=True,
  93.                                marker=dict(size=10, opacity=1, color='black', symbol=donor_symbol),
  94.                                hoverinfo='none'))
  96. # Annotations for the legend blocks
  97. fig.update_layout(
  98.     annotations=[
  99.         dict(x=1.1, y=1.0, xref='paper', yref='paper', showarrow=False,
  100.              text='Condition', font=dict(size=15)),
  101.         dict(x=1.1, y=0.6, xref='paper', yref='paper', showarrow=False,
  102.              text='Donor', font=dict(size=15))
  103.     ],
  104.     scene=dict(
  105.         aspectmode='cube',
  106.         xaxis=dict(gridcolor='black', backgroundcolor='white', zerolinecolor='black', title='PC1: 36% v.'),
  107.         yaxis=dict(gridcolor='black', backgroundcolor='white', zerolinecolor='black', title='PC2: 17% v.'),
  108.         zaxis=dict(gridcolor='black', backgroundcolor='white', zerolinecolor='black', title='PC3: 15% variance'),
  109.         bgcolor='white'
  110.     ),
  111.     margin=dict(l=5, r=5, b=5, t=0)  # Adjust the margins to prevent clipping of axis titles
  112. )
  114. #fig.show()
  115. fig.write_image("fig1.svg")

calculate background genes

  1.     # [1] ensembl_gene_id         gene_name               MKL.1_RNA              
  2.     # [4] MKL.1_RNA_118           MKL.1_RNA_147           MKL.1_EV.RNA           
  3.     # [7] MKL.1_EV.RNA_2          MKL.1_EV.RNA_118        MKL.1_EV.RNA_87        
  4.     #[10] MKL.1_EV.RNA_27         X042_MKL.1_wt_EV        X0505_MKL.1_wt_EV      
  5.     #[13] X042_MKL.1_sT_DMSO      X0505_MKL.1_sT_DMSO_EV  X042_MKL.1_scr_DMSO_EV 
  6.     #[16] X0505_MKL.1_scr_DMSO_EV X042_MKL.1_sT_Dox       X0505_MKL.1_sT_Dox_EV  
  7.     #[19] X042_MKL.1_scr_Dox_EV   X0505_MKL.1_scr_Dox_EV
  9. "ensembl_gene_id
  10. gene_name 
  11. python process.py ../fpkm_values_MKL-1.csv 3 > MKL.1_RNA 
  12. python process.py ../fpkm_values_MKL-1.csv 4 > MKL.1_RNA_118 
  13. python process.py ../fpkm_values_MKL-1.csv 5 > MKL.1_RNA_147 
  14. python process.py ../fpkm_values_MKL-1.csv 6 > MKL.1_EV.RNA 
  15. python process.py ../fpkm_values_MKL-1.csv 7 > MKL.1_EV.RNA_2 
  16. python process.py ../fpkm_values_MKL-1.csv 8 > MKL.1_EV.RNA_118 
  17. python process.py ../fpkm_values_MKL-1.csv 9 > MKL.1_EV.RNA_87 
  18. python process.py ../fpkm_values_MKL-1.csv 10 > MKL.1_EV.RNA_27 
  19. python process.py ../fpkm_values_MKL-1.csv 11 > X042_MKL.1_wt_EV 
  20. python process.py ../fpkm_values_MKL-1.csv 12 > X0505_MKL.1_wt_EV 
  21. python process.py ../fpkm_values_MKL-1.csv 13 > X042_MKL.1_sT_DMSO 
  22. python process.py ../fpkm_values_MKL-1.csv 14 > X0505_MKL.1_sT_DMSO_EV 
  23. python process.py ../fpkm_values_MKL-1.csv 15 > X042_MKL.1_scr_DMSO_EV 
  24. python process.py ../fpkm_values_MKL-1.csv 16 > X0505_MKL.1_scr_DMSO_EV 
  25. python process.py ../fpkm_values_MKL-1.csv 17 > X042_MKL.1_sT_Dox 
  26. python process.py ../fpkm_values_MKL-1.csv 18 > X0505_MKL.1_sT_Dox_EV 
  27. python process.py ../fpkm_values_MKL-1.csv 19 > X042_MKL.1_scr_Dox_EV 
  28. python process.py ../fpkm_values_MKL-1.csv 20 > X0505_MKL.1_scr_Dox_EV 
  29. python process.py ../fpkm_values_MKL-1.csv 21 > median_length
  31. cut -f1 -d',' MKL.1_RNA>_MKL-1_RNA
  32. cut -f1 -d',' MKL.1_RNA_118>_MKL-1_RNA_118
  33. cut -f1 -d',' MKL.1_RNA_147>_MKL-1_RNA_147
  34. cut -f1 -d',' MKL.1_EV.RNA>_MKL-1_EV_RNA
  35. cut -f1 -d',' MKL.1_EV.RNA_2>_MKL-1_EV_RNA_2
  36. cut -f1 -d',' MKL.1_EV.RNA_118>_MKL-1_EV_RNA_118
  37. cut -f1 -d',' MKL.1_EV.RNA_87>_MKL-1_EV_RNA_87
  38. cut -f1 -d',' MKL.1_EV.RNA_27>_MKL-1_EV_RNA_27
  39. cut -f1 -d',' X042_MKL.1_wt_EV>_X042_MKL-1_wt_EV
  40. cut -f1 -d',' X0505_MKL.1_wt_EV>_X0505_MKL-1_wt_EV
  41. cut -f1 -d',' X042_MKL.1_sT_DMSO>_X042_MKL-1_sT_DMSO
  42. cut -f1 -d',' X0505_MKL.1_sT_DMSO_EV>_X0505_MKL-1_sT_DMSO_EV
  43. cut -f1 -d',' X042_MKL.1_scr_DMSO_EV>_X042_MKL-1_scr_DMSO_EV
  44. cut -f1 -d',' X0505_MKL.1_scr_DMSO_EV>_X0505_MKL-1_scr_DMSO_EV
  45. cut -f1 -d',' X042_MKL.1_sT_Dox>_X042_MKL-1_sT_Dox
  46. cut -f1 -d',' X0505_MKL.1_sT_Dox_EV>_X0505_MKL-1_sT_Dox_EV
  47. cut -f1 -d',' X042_MKL.1_scr_Dox_EV>_X042_MKL-1_scr_Dox_EV
  48. cut -f1 -d',' X0505_MKL.1_scr_Dox_EV>_X0505_MKL-1_scr_Dox_EV
  50. ~/Tools/csv2xls-0.4/csv_to_xls.py  _MKL-1_RNA _MKL-1_RNA_118 _MKL-1_RNA_147 _MKL-1_EV_RNA _MKL-1_EV_RNA_2 _MKL-1_EV_RNA_118 _MKL-1_EV_RNA_87 _MKL-1_EV_RNA_27 _X042_MKL-1_wt_EV _X0505_MKL-1_wt_EV _X042_MKL-1_sT_DMSO _X0505_MKL-1_sT_DMSO_EV _X042_MKL-1_scr_DMSO_EV _X0505_MKL-1_scr_DMSO_EV _X042_MKL-1_sT_Dox _X0505_MKL-1_sT_Dox_EV _X042_MKL-1_scr_Dox_EV _X0505_MKL-1_scr_Dox_EV -d',' -o background_genes_MKL-1.xls;

process.py

  1. import sys
  3. filename = sys.argv[1]
  4. n = int(sys.argv[2])
  6. with open(filename, 'r') as file:
  7.     for line in file:
  8.         words = line.split(",")
  9.         if len(words) >= n and words[n-1] != '0':
  10.             print(line.strip())
  12. # -- in the example of K331A data --
  13. #FORMAT normalized_counts.txt to normalized_counts.csv
  14. #DEBUG1: add gene_symbol in the first line
  15. #DEBUG2: '\n' --> ',\n'
  16. #RESULT: it looks like #gene_symbol,ctrl_d3_DonorI,ctrl_d3_DonorII,ctrl_d8_DonorI,ctrl_d8_DonorII,LT_d3_DonorI,LT_d3_DonorII,LT_d8_DonorI,LT_d8_DonorII,LTtr_d3_DonorI,LTtr_d3_DonorII,LTtr_d8_DonorI,LTtr_d8_DonorII,K331A_DonorI,K331A_DonorII,
  17. #DDX11L1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  18. python process.py normalized_counts.csv 2 > ctrl_d3_DonorI
  19. python process.py normalized_counts.csv 3 > ctrl_d3_DonorII 
  20. python process.py normalized_counts.csv 4 > ctrl_d8_DonorI
  21. python process.py normalized_counts.csv 5 > ctrl_d8_DonorII
  22. python process.py normalized_counts.csv 6 > LT_d3_DonorI
  23. python process.py normalized_counts.csv 7 > LT_d3_DonorII
  24. python process.py normalized_counts.csv 8 > LT_d8_DonorI
  25. python process.py normalized_counts.csv 9 > LT_d8_DonorII
  26. python process.py normalized_counts.csv 10 > LTtr_d3_DonorI
  27. python process.py normalized_counts.csv 11 > LTtr_d3_DonorII
  28. python process.py normalized_counts.csv 12 > LTtr_d8_DonorI
  29. python process.py normalized_counts.csv 13 > LTtr_d8_DonorII
  30. python process.py normalized_counts.csv 14 > K331A_DonorI
  31. python process.py normalized_counts.csv 15 > K331A_DonorII
  33. cut -f1 -d',' ctrl_d3_DonorI > _ctrl_d3_DonorI
  34. cut -f1 -d',' ctrl_d3_DonorII > _ctrl_d3_DonorII
  35. cut -f1 -d',' ctrl_d8_DonorI > _ctrl_d8_DonorI
  36. cut -f1 -d',' ctrl_d8_DonorII > _ctrl_d8_DonorII
  37. cut -f1 -d',' LT_d3_DonorI > _LT_d3_DonorI
  38. cut -f1 -d',' LT_d3_DonorII > _LT_d3_DonorII
  39. cut -f1 -d',' LT_d8_DonorI > _LT_d8_DonorI
  40. cut -f1 -d',' LT_d8_DonorII > _LT_d8_DonorII
  41. cut -f1 -d',' LTtr_d3_DonorI > _LTtr_d3_DonorI
  42. cut -f1 -d',' LTtr_d3_DonorII > _LTtr_d3_DonorII
  43. cut -f1 -d',' LTtr_d8_DonorI > _LTtr_d8_DonorI
  44. cut -f1 -d',' LTtr_d8_DonorII > _LTtr_d8_DonorII
  45. cut -f1 -d',' K331A_DonorI > _K331A_DonorI
  46. cut -f1 -d',' K331A_DonorII > _K331A_DonorII
  48. ~/Tools/csv2xls-0.4/csv_to_xls.py _ctrl_d3_DonorI _ctrl_d3_DonorII _ctrl_d8_DonorI _ctrl_d8_DonorII   _LT_d3_DonorI _LT_d3_DonorII _LT_d8_DonorI _LT_d8_DonorII  _LTtr_d3_DonorI _LTtr_d3_DonorII _LTtr_d8_DonorI _LTtr_d8_DonorII  _K331A_DonorI _K331A_DonorII  -d',' -o background_genes.xls;

like unlike

点赞本文的读者

还没有人对此文章表态

本文有评论

没有评论

看文章，发评论，不要沉默

最受欢迎文章

  1. Why Do Significant Gene Lists Change After Adding Additional Conditions in Differential Gene Expression Analysis?
  2. Motif Discovery in Biological Sequences: A Comparison of MEME and HOMER
  3. Updating Human Gene Identifiers using Ensembl BioMart: A Step-by-Step Guide
  4. Setup conda environments
  5. Calling peaks using findPeaks of HOMER
  6. File format for single channel analysis of Agilent microarray data with Limma?
  7. pheatmap vs heatmap.2
  8. Should the inputs for GSVA be normalized or raw?
  9. 洛那法尼(lonafarnib)：从抗癌症到抗病毒的多功能药物
  10. RNAseq running with umi_tools

最新文章

最多评论文章

推荐相似文章

YopQ Secretion Boxplot and Fitting Function

RNA-seq skin organoids on GRCh38+chrHsv1 (final)

Genomic Organization of Herpes Simplex Virus type 1 (HSV-1 s17)

From epidome to treeHeatmap

© 2023 XGenes.com Impressum