{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"Here I checked the overlap ration between CITE-seq input genes and ATAC-seq target genes","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"!pip install -q tables","metadata":{"execution":{"iopub.status.busy":"2022-09-01T02:02:19.890610Z","iopub.execute_input":"2022-09-01T02:02:19.891180Z","iopub.status.idle":"2022-09-01T02:02:33.823922Z","shell.execute_reply.started":"2022-09-01T02:02:19.891078Z","shell.execute_reply":"2022-09-01T02:02:33.822466Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n%matplotlib inline\nimport matplotlib.pyplot as plt\nfrom matplotlib_venn import venn2\nimport seaborn as sns\nsns.set()","metadata":{"execution":{"iopub.status.busy":"2022-09-01T02:19:56.724871Z","iopub.execute_input":"2022-09-01T02:19:56.725528Z","iopub.status.idle":"2022-09-01T02:19:57.077513Z","shell.execute_reply.started":"2022-09-01T02:19:56.725483Z","shell.execute_reply":"2022-09-01T02:19:57.076285Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cite_data = pd.read_hdf('../input/open-problems-multimodal/train_cite_inputs.h5', start=0, stop=5)\ncite_data.head()","metadata":{"execution":{"iopub.status.busy":"2022-09-01T02:02:34.034768Z","iopub.execute_input":"2022-09-01T02:02:34.035225Z","iopub.status.idle":"2022-09-01T02:02:34.223493Z","shell.execute_reply.started":"2022-09-01T02:02:34.035180Z","shell.execute_reply":"2022-09-01T02:02:34.222316Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cite_genes = {gene.split('_')[0] for gene in cite_data.columns}","metadata":{"execution":{"iopub.status.busy":"2022-09-01T02:02:34.226779Z","iopub.execute_input":"2022-09-01T02:02:34.227915Z","iopub.status.idle":"2022-09-01T02:02:34.245758Z","shell.execute_reply.started":"2022-09-01T02:02:34.227867Z","shell.execute_reply":"2022-09-01T02:02:34.244810Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"atac_data = pd.read_hdf('../input/open-problems-multimodal/train_multi_targets.h5', start=0, stop=5)\natac_data.head()","metadata":{"execution":{"iopub.status.busy":"2022-09-01T02:02:34.247313Z","iopub.execute_input":"2022-09-01T02:02:34.248080Z","iopub.status.idle":"2022-09-01T02:02:34.382329Z","shell.execute_reply.started":"2022-09-01T02:02:34.248036Z","shell.execute_reply":"2022-09-01T02:02:34.381178Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"atac_genes = set(atac_data.columns)","metadata":{"execution":{"iopub.status.busy":"2022-09-01T02:02:34.383586Z","iopub.execute_input":"2022-09-01T02:02:34.383895Z","iopub.status.idle":"2022-09-01T02:02:34.394140Z","shell.execute_reply.started":"2022-09-01T02:02:34.383867Z","shell.execute_reply":"2022-09-01T02:02:34.392627Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"diff_cite = cite_genes - atac_genes\ndiff_atac = atac_genes - cite_genes\ninter = cite_genes & atac_genes\nvenn2(subsets = (len(diff_cite), len(diff_atac), len(inter)), set_labels = ('Genes in CITE inputs', 'Genes in ATAC targets'))\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-01T02:20:00.951053Z","iopub.execute_input":"2022-09-01T02:20:00.951460Z","iopub.status.idle":"2022-09-01T02:20:01.075531Z","shell.execute_reply.started":"2022-09-01T02:20:00.951427Z","shell.execute_reply":"2022-09-01T02:20:01.074326Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"CITE-seq genes and ATAC-seq genes had considerable overlap(around 80%).","metadata":{}},{"cell_type":"markdown","source":"## RNA annotations\nRef: https://rnact.crg.eu/","metadata":{}},{"cell_type":"code","source":"rna_annot = pd.read_table('../input/rnact-supporting-tables/catrapid_rnas.txt')\nrna_annot.head()","metadata":{"execution":{"iopub.status.busy":"2022-09-01T02:03:46.603189Z","iopub.execute_input":"2022-09-01T02:03:46.603688Z","iopub.status.idle":"2022-09-01T02:03:48.575706Z","shell.execute_reply.started":"2022-09-01T02:03:46.603644Z","shell.execute_reply":"2022-09-01T02:03:48.574595Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rna_annot_human = rna_annot[rna_annot['species'] == 'human'].reset_index(drop=True)\ndisplay(rna_annot_human.shape)","metadata":{"execution":{"iopub.status.busy":"2022-09-01T02:04:47.061574Z","iopub.execute_input":"2022-09-01T02:04:47.062250Z","iopub.status.idle":"2022-09-01T02:04:47.218077Z","shell.execute_reply.started":"2022-09-01T02:04:47.062213Z","shell.execute_reply":"2022-09-01T02:04:47.217184Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rna_biotype_dict = {ens: biotype for ens, biotype in zip(rna_annot_human['ensg'], rna_annot_human['biotype'])}","metadata":{"execution":{"iopub.status.busy":"2022-09-01T02:15:45.854230Z","iopub.execute_input":"2022-09-01T02:15:45.854668Z","iopub.status.idle":"2022-09-01T02:15:45.927962Z","shell.execute_reply.started":"2022-09-01T02:15:45.854633Z","shell.execute_reply":"2022-09-01T02:15:45.926827Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Overlap genes","metadata":{}},{"cell_type":"code","source":"# genes not registered in RNAct database\nexcl_genes = inter - set(rna_annot_human['ensg'])\nprint(f'{len(excl_genes)} genes ({len(excl_genes)/len(inter)*100:.2f}%) are not registered in RNAct database.')\nfor n, gene in enumerate(excl_genes):\n    print(gene.ljust(16), end='')\n    if (n+1) % 10 == 0: print('\\n')","metadata":{"execution":{"iopub.status.busy":"2022-09-01T02:16:54.763472Z","iopub.execute_input":"2022-09-01T02:16:54.763859Z","iopub.status.idle":"2022-09-01T02:16:54.802993Z","shell.execute_reply.started":"2022-09-01T02:16:54.763829Z","shell.execute_reply":"2022-09-01T02:16:54.801638Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"incl_genes = inter & set(rna_annot_human['ensg'])\noverlap_genes_biotype = [rna_biotype_dict[ens] for ens in incl_genes]\nfig = plt.figure(figsize=(6, 7))\nax = fig.add_subplot()\npd.Series(overlap_genes_biotype).value_counts(ascending=True).plot.barh(ax=ax)\nax.set_title('Genes found both in CITE-seq inputs and in ATAC-seq targets\\n')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-01T02:34:16.154469Z","iopub.execute_input":"2022-09-01T02:34:16.154886Z","iopub.status.idle":"2022-09-01T02:34:16.571679Z","shell.execute_reply.started":"2022-09-01T02:34:16.154841Z","shell.execute_reply":"2022-09-01T02:34:16.570431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Most of overlap genes were categorized as protein coding.","metadata":{}},{"cell_type":"markdown","source":"### Genes found only in CITE-seq inputs","metadata":{"execution":{"iopub.status.busy":"2022-09-01T02:31:36.690298Z","iopub.execute_input":"2022-09-01T02:31:36.690794Z","iopub.status.idle":"2022-09-01T02:31:37.105003Z","shell.execute_reply.started":"2022-09-01T02:31:36.690757Z","shell.execute_reply":"2022-09-01T02:31:37.103690Z"}}},{"cell_type":"code","source":"# genes not registered in RNAct database\nexcl_genes = diff_cite - set(rna_annot_human['ensg'])\nprint(f'{len(excl_genes)} genes ({len(excl_genes)/len(diff_cite)*100:.2f}%) are not registered in RNAct database.')\nfor n, gene in enumerate(excl_genes):\n    print(gene.ljust(16), end='')\n    if (n+1) % 10 == 0: print('\\n')","metadata":{"execution":{"iopub.status.busy":"2022-09-01T02:36:20.287096Z","iopub.execute_input":"2022-09-01T02:36:20.288223Z","iopub.status.idle":"2022-09-01T02:36:20.324757Z","shell.execute_reply.started":"2022-09-01T02:36:20.288175Z","shell.execute_reply":"2022-09-01T02:36:20.323494Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"incl_genes = diff_cite & set(rna_annot_human['ensg'])\noverlap_genes_biotype = [rna_biotype_dict[ens] for ens in incl_genes]\nfig = plt.figure(figsize=(6, 7))\nax = fig.add_subplot()\npd.Series(overlap_genes_biotype).value_counts(ascending=True).plot.barh(ax=ax)\nax.set_title('Genes found only in CITE-seq inputs\\n')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-01T02:34:51.793282Z","iopub.execute_input":"2022-09-01T02:34:51.793705Z","iopub.status.idle":"2022-09-01T02:34:52.265146Z","shell.execute_reply.started":"2022-09-01T02:34:51.793674Z","shell.execute_reply":"2022-09-01T02:34:52.264066Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Most of genes found only in CITE-seq inputs were categorized as pseudogene.","metadata":{}},{"cell_type":"markdown","source":"### Genes found only in ATAC-seq targets","metadata":{}},{"cell_type":"code","source":"# genes not registered in RNAct database\nexcl_genes = diff_atac - set(rna_annot_human['ensg'])\nprint(f'{len(excl_genes)} genes ({len(excl_genes)/len(diff_atac)*100:.2f}%) are not registered in RNAct database.')\nfor n, gene in enumerate(excl_genes):\n    print(gene.ljust(16), end='')\n    if (n+1) % 10 == 0: print('\\n')","metadata":{"execution":{"iopub.status.busy":"2022-09-01T02:36:55.984939Z","iopub.execute_input":"2022-09-01T02:36:55.985553Z","iopub.status.idle":"2022-09-01T02:36:56.038094Z","shell.execute_reply.started":"2022-09-01T02:36:55.985511Z","shell.execute_reply":"2022-09-01T02:36:56.037027Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"18.3% of genes found only in ATAC-seq targets were not well known and it was higher than CITE-seq case.","metadata":{}},{"cell_type":"code","source":"incl_genes = diff_atac & set(rna_annot_human['ensg'])\noverlap_genes_biotype = [rna_biotype_dict[ens] for ens in incl_genes]\nfig = plt.figure(figsize=(6, 8))\nax = fig.add_subplot()\npd.Series(overlap_genes_biotype).value_counts(ascending=True).plot.barh(ax=ax)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-01T02:37:07.652888Z","iopub.execute_input":"2022-09-01T02:37:07.653367Z","iopub.status.idle":"2022-09-01T02:37:08.077546Z","shell.execute_reply.started":"2022-09-01T02:37:07.653306Z","shell.execute_reply":"2022-09-01T02:37:08.076218Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Most of genes found only in ATAC-seq inputs were categorized as long noncoding RNA including antisense RNA.","metadata":{}},{"cell_type":"markdown","source":"* CITE-seq input genes and ATAC-seq target genes had considerable overlap (around 80%).\n* Most of overlap genes were categorized as protein coding.\n* Most of genes found only in CITE-seq inputs were categorized as pseudogene. (Should we drop them?)\n* 18.3% of genes found only in ATAC-seq targets were not well known and it was higher than CITE-seq case.\n* Most of genes found only in ATAC-seq inputs were categorized as long noncoding RNA including antisense RNA.\n* (The prediction of noncoding RNA expression might be difficult.)","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}