{"cells":[{"metadata":{},"cell_type":"markdown","source":"#
Deeper Exploratory Data Analysis of PANDA
\n\n\n## About the notebook\n\nIn this notebook, we are basically going to dive deeper into the dataset explanation and build on that to perform an extensive EDA. However, to get started we will first begin with a deep explanation of some key concepts related to **Prostate Cancer** and its detection and that we are worth knowing.\n","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Domain Knowledge \n\nIn this section, we are interested in covering the knowledge around the Prostate Cancer in order to have a deeper understanding of the problem. To do this, we are going to start by trying to successively address these five key questions: what is Prostate Cancer? How is it tested and detected? What is Gleason score and how does it fit into all this? What is ISUP grade and how is it related to Gleason score? How have Gleason scores been generated?\n\nSo, let us start looking into these questions\n\n### Question 1: What is Prostate Cancer?\n\nProstate Cancer is a cancer that occurs in the prostate. The prostate is nothing but a small walnut-shaped gland in men that produces the seminal fluid that nourishes and transports sperm as depicted in the figure below.\n\n
\n\n \n \n \n \n
\n \n
\n
\n Prostate Cancer is one of the most common types of cancer in men. Usually, it grows slowly and is initially confined to the prostate gland, where it may not cause serious harm. However, while some types of Prostate Cancer grow slowly and may need minimal or even no treatment, other types are aggressive and can spread quickly.\n
\n
\n
\n\n\n","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Question 2: How is Prostate Cancer tested and detected?\n\nProstate Cancer is tested through some prostate screening tests that might include:\n* **Digital rectal exam (DRE)**: During a DRE, your doctor inserts a gloved, lubricated finger into your rectum During a DRE, your doctor inserts a gloved, lubricated finger into your rectum to examine your prostate, which is adjacent to the rectum as depicted in the figure below. More precisely, he feels the back wall of the prostate gland for enlargement, tenderness, lumps or hard spots. If your doctor finds any abnormalities in the texture, shape or size of the gland, you may need further tests.\n
\n\n \n \n \n
\n
\n* **Prostate-specific antigen (PSA) test**: Here, a blood sample is drawn from a vein in your arm and analyzed for PSA, a substance that's naturally produced by your prostate gland. It's normal for a small amount of PSA to be in your bloodstream. However, if a higher than normal level is found, it may indicate prostate infection, inflammation, enlargement or cancer.\n\nIf a DRE or PSA test detects an abnormality, your doctor may recommend further tests to determine whether you have prostate cancer. Such tests can include:\n\n* **Ultrasound**: If other tests raise concerns, your doctor may use transrectal ultrasound to further evaluate your prostate. A small probe, about the size and shape of a cigar, is inserted into your rectum. The probe uses sound waves to create a picture of your prostate gland.\n* **Collecting a sample of prostate tissue** : If initial test results suggest prostate cancer, your doctor may recommend a procedure to collect a sample of cells from your prostate (prostate biopsy). Prostate biopsy is often done using a thin needle that's inserted into the prostate to collect tissue. The tissue sample is analyzed in a lab to determine whether cancer cells are present.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Question 3: What is Gleason score and how does it fit into all this?\n\nWhen a biopsy confirms the presence of cancer, the next step is to determine the level of aggressiveness (grade) of the cancer cells. A laboratory pathologist examines a sample of your cancer to determine how much cancer cells differ from the healthy cells. A higher grade indicates a more aggressive cancer that is more likely to spread quickly.\n\nThe most common scale used to evaluate the grade of prostate cancer cells is called a Gleason score. Gleason scoring combines two numbers and can range from 1-10 and describes how much the cancer from a biopsy looks like healthy tissue (lower score or less aggressive) or abnormal tissue (higher score or very aggressive)., though the lower part of the range isn't used as often.\n\n### Question 4: What is ISUP grade and how is it related to Gleason score?\n\nAccording to current guidelines by the International Society of Urological Pathology (ISUP), the Gleason scores are summarized into an ISUP grade on a scale from 1 to 5 according to the following rule:\n* Gleason score 6 = ISUP grade 1\n* Gleason score 7 (3 + 4) = ISUP grade 2\n* Gleason score 7 (4 + 3) = ISUP grade 3\n* Gleason score 8 = ISUP grade 4\n* Gleason score 9-10 = ISUP grade 5\n\nIf there is no cancer in the sample, we use the label ISUP grade 0 in this competition. ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Question 5: How have Gleason scores been generated in the dataset?\n\nEach whole-slide Image (WSI) in this challenge contains one, or in some cases two, thin tissue sections cut from a single biopsy sample. Prior to scanning, the tissue is stained with haematoxylin & eosin (H&E). This is a standard way of staining the originally transparent tissue to produce some contrast. The samples are made up of glandular tissue and connective tissue. The glands are hollow structures, which can be seen as white “holes” or branched cavities in the WSI. The appearance of the glands forms the basis of the Gleason grading system. The glandular structure characteristic of healthy prostate tissue is progressively lost with increasing grade. The grading system recognizes three categories: 3, 4, and 5. The patterns are described in detail below and exemplified in the figure below:\n\n* **[A] Benign prostate glands with folded epithelium:** The cytoplasm is pale and the nuclei small and regular. The glands are grouped together.\n* **[B] Prostatic adenocarcinoma:** Gleason Pattern 3 has no loss of glandular differentiation. Small glands infiltrate between benign glands. The cytoplasm is often dark and the nuclei enlarged with dark chromatin and some prominent nucleoli. Each epithelial unit is separate and has a lumen.\n* **[C] Prostatic adenocarcinoma:** Gleason Pattern 4 has partial loss of glandular differentiation. There is an attempt to form lumina but the tumor fails to form complete, well-developed glands. This microphotograph shows irregular cribriform cancer, i.e. epithelial sheets with multiple lumina. There are also some poorly formed small glands and some fused glands. All of these are included in Gleason Pattern 4.\n* **[D] Prostatic adenocarcinoma:** Gleason Pattern 5 has an almost complete loss of glandular differentiation. Dispersed single cancer cells are seen in the stroma. Gleason Pattern 5 may also contain solid sheets or strands of cancer cells. All microphotographs show hematoxylin and eosin stains at 20x lens magnification.\n\n
\n\n \n \n \n \n
\n
\n Examples of patterns used for grading. [A] Healthy glands. [B]-[D] Gleason pattern 3 to 5, respectively.\n
\n
\n
\n","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Now Let's Get Started with PANDA Dataset\n\nThis notebook shows few methods to load and preprocess WSIs from PANDA challenge dataset. Base on the challenge description, the dataset consists of about 11 000 data samples for which each sample represents a whole-slide image (WSI) of prostate biopsies from two sources: [Radboud University Medical Center](https://www.radboudumc.nl/en/research) and the [Karolinska Institute](https://ki.se/en/meb).\n\n\n## Loading of the Required Dependencies\n\nThis Python 3 environment associated with this notebook comes with many helpful analytics libraries installed. It is defined by the `kaggle/python` Docker image: [https://github.com/kaggle/docker-python]() In the cell below we enable the several different and helpful packages/modules/libraries that I think are required for this project.","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"### Directory dependencies\nimport os\nfrom glob import glob\n\n### Preprocessing dependencies: Level 1\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport random\n\n### Plotting dependencies\nimport matplotlib as mpl\nimport matplotlib.pyplot as plt\nimport plotly.graph_objs as go # Plotly for the interactive viewer (see last section)\nimport seaborn as sns\nfrom IPython.display import Image, display, HTML\n\n### Preprocessing dependencies: Level 2\nimport openslide as OS\nimport PIL\nfrom PIL import Image, ImageOps\nimport imageio\nimport cv2\nimport skimage.io as IO\n\n### Figure Configuration\nsns.set_style(\"whitegrid\")\nplt.rc(\"figure\", titlesize=20) # fontsize of the figure title\nplt.rc(\"axes\", titlesize=17) # fontsize of the axes title\nplt.rc(\"axes\", labelsize=15) # fontsize of the x and y labels\nplt.rc(\"xtick\", labelsize=12) # fontsize of the tick labels\nplt.rc(\"ytick\", labelsize=12) # fontsize of the tick labels\nplt.rc(\"legend\", fontsize=13) # legend fontsize","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Loading of the Data\n\n### Directories' Exploration\n\n**Note:** We can write up to 5GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\". You can also write temporary files to `/kaggle/temp/`, but they won't be saved outside of the current session.","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"pwd","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nprint(\"This is an overview of the different directories and their contents\\n\")\ncount = 1\nfor dirname, _, filenames in os.walk(\"/kaggle/\"):\n for filename in filenames:\n print(os.path.join(dirname, filename))\n count += 1\n \n if count%10 == 0:\n break","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Load the Data and Check the Content","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"### Base directory of the data\ndata_dir = \"../input/prostate-cancer-grade-assessment\"\n\n### Directories of the training images and label masks\ntrain_dir = os.path.sep.join([data_dir, \"train_images\"])\ntrain_mask_dir = os.path.sep.join([data_dir, \"train_label_masks\"])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"### Load the data\ntrain_ids = pd.read_csv(f\"{data_dir}/train.csv\")\ntest_ids = pd.read_csv(f\"{data_dir}/test.csv\")\nsubmission = pd.read_csv(f\"{data_dir}/sample_submission.csv\")\n\ndisplay(train_ids.tail(3).style.background_gradient(cmap=\"Blues\"))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"train_ids.set_index(\"image_id\", inplace=True)\ndisplay(train_ids.tail(10).style.background_gradient(cmap=\"Blues\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Checking for the Training Samples","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"###\nfiles = glob(f\"{train_dir}/*\")\nprint(f\"Size of the samples in the train_images directory: {len(files)}\")\nprint(f\"Size of training samples whose IDs are in train.csv: {train_ids.shape[0]}\")\nfiles = glob(f\"{train_mask_dir}/*\")\nprint(f\"Size of the samples masks in the train_label_masks directory: {len(files)}\")\n###\nprint(f\"\\nData providers IDs: {train_ids.data_provider.unique()}, Size: {len(train_ids.data_provider.unique())}\")\nprint(f\"ISUP Grades (target) used: {train_ids.isup_grade.unique()}, Size: {len(train_ids.isup_grade.unique())}\")\nprint(f\"Gleason Score used: {train_ids.gleason_score.unique()}, Size: {len(train_ids.gleason_score.unique())}\")\n###\n#display(train_ids.tail(10).style.background_gradient(cmap=\"Blues\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Inference\n\nWe can notice from the outputs of the cell above that\n* The training samples directory, `train_images/`, contains exactly 10616 files that matches the number of samples in the `train.csv`.\n* All the training images do not have masks, as we can see that the size of the training samples is **10616** whereas the sample mask size is only **10516**. This means there are exactly 100 training samples without that do not have a mask.\n\n#### About the Output Variables\n\nThe whole experiment of testing and detecting the Prostate Cancer outputs two varaibles that are:\n\n* **isup_grade**: which is the target variable representing the severity of the cancer on a `0-5` scale.\n* **gleason_score**: which is an alternate cancer severity rating system with more levels than the ISUP scale as you can see in the figure below depicting how Gleason score and ISUP grade systems compare. Let's recall this is provided on training data only.\n\n
\n\n \n \n \n \n
\n
\n An illustration of the Gleason grading process for an example biopsy containing prostate cancer. The most common (blue outline, Gleason pattern 3) and second most common (red outline, Gleason pattern 4) cancer growth patterns present in the biopsy dictate the Gleason score (3+4 for this biopsy), which in turn is converted into an ISUP grade (2 for this biopsy) following guidelines of the International Society of Urological Pathology. Biopsies not containing cancer are represented by an ISUP grade of 0 in this challenge.\n
\n
\n
\n\n\n#### Checking Missing Information","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"train_ids.isnull().sum()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Checking of `test.csv` File Content","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"display(test_ids.head())\n###\nprint(f\"Size of test samples whose IDs are in test.csv: {test_ids.shape[0]}\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"# Performing of Some Basic EDA\n\n## Checking of the Size of `Training label masks` and Samples and \n\n","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"df = pd.DataFrame({\"data_type\": [\"train_images\", \"train_label_masks\"],\n \"count\": [len(glob(f\"{train_dir}/*\")), len(glob(f\"{train_mask_dir}/*\"))]}).set_index(\"data_type\")\n\ndisplay(df.style.background_gradient(cmap=\"Blues\"))\n####\nfig, ax = plt.subplots(figsize=[6,8])\n\nval_ax = sns.barplot(x=df.index, y=df[\"count\"], palette=\"deep\", ax=ax)\nfor i, v in enumerate(df.values):\n ax.text(val_ax.get_xticks()[i], v, str(int(v)),\n ha=\"center\", fontsize=13)\nax.set_ylabel(\"count\")\nax.set_title(\"Training sample Size\\n\")\nax.set_xlabel(\"\\ndata type\")\nfig.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Inference\n\nAs we also highlighted early in the previous section, the outputs of the cell above show that:\n* All the training images do not have masks, as we can see that the training samples size is **10616** whereas the sample mask size is only **10516**.\n\n## Checking of the Distribution of `Data Provider`, `ISUP Grade` and `Gleason Score`\n\n### Auxilary Function","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"def distribution_plot(df, feature, ax, title=\"\"):\n total = float(len(df))\n sns.countplot(df[feature],\n order=df[feature].value_counts().index,\n ax=ax)\n \n if feature==\"gleason_score\":\n ax.set_xticklabels(df[feature].tolist(), rotation=45) \n ax.set_title(title)\n if feature!=\"data_provider\":\n ax.set_ylabel(\"\")\n for p in ax.patches:\n height = p.get_height()\n ax.text(p.get_x()+p.get_width()/2.,\n height+3,\n \"{:1.1f}%\".format(100*height/total),\n ha=\"center\",\n fontsize=12) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Summary of Counts","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"df = train_ids.reset_index()\ntemp1 = df.groupby(\"data_provider\").count()[\"image_id\"].reset_index().sort_values(by=\"image_id\",ascending=False).rename(columns={\"image_id\": \"count\"})\ntemp2 = df.groupby(\"isup_grade\").count()[\"image_id\"].reset_index().sort_values(by=\"image_id\", ascending=False).rename(columns={\"image_id\": \"count\"})\ntemp3 = df.groupby(\"gleason_score\").count()[\"image_id\"].reset_index().sort_values(by=\"image_id\",ascending=False).rename(columns={\"image_id\": \"count\"})\ndisplay(temp1.style.background_gradient(cmap=\"Blues\"))\ndisplay(temp2.style.background_gradient(cmap=\"Blues\"))\ndisplay(temp3.style.background_gradient(cmap=\"Blues\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Visualization","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"df = train_ids.reset_index()\nfig = plt.figure(figsize=(22, 8))\nax = [fig.add_subplot(1, 3, 1),\n fig.add_subplot(1, 3, 2),\n fig.add_subplot(1, 3, 3)]\n###\ntitle0 = \"Data provider\"\ntitle1 = \"ISUP grade\"\ntitle2 = \"Gleason score\"\n###\ndistribution_plot(df=df, feature=\"data_provider\", ax=ax[0], title=title0)\ndistribution_plot(df=df, feature=\"isup_grade\", ax=ax[1], title=title1)\ndistribution_plot(df=df, feature=\"gleason_score\", ax=ax[2], title=title2)\n###\nfig.suptitle(\"Distribution plot (count and %)\", y=1.1)\nfig.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Inferences on the Distributions\n\n#### About the Distribution of `ISUP Grade`\n\n* Majority of data samples in train set have ISUP grade values 0 or 1 (total > 50%).\n* Rest of the data samples have associated ISUP grades from 2 to 5 with all ranging in the 11-12% each.\n* Dataset is not balanced in terms of isup_grade.\n\n\n#### About the Distribution of `Gleason Socre`\n\n* From above graph, it is clear that gleason_score distribution is not uniform.\n* Few gleason_score like (3+3) and (0+0) are more frequent while others like (3+5) and (5+3) are very rare in this dataset.\n* Dataset is not balanced in terms of gleason_score.\n\n\n## Checking for the Relative Distributions (RDs)\n1. `ISUP Grade` and `Data Provider`\n1. `Gleason Score` and `Data Provider`\n","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"df = train_ids.reset_index()\nfig, ax = plt.subplots(1, 2, figsize=[20,8], sharey=True)\n\nval1 = df.groupby([\"isup_grade\", \"data_provider\"]).count()[\"image_id\"].unstack().plot(kind=\"bar\", ax=ax[0])\nval2 = df.groupby([\"gleason_score\", \"data_provider\"]).count()[\"image_id\"].unstack().plot(kind=\"bar\", ax=ax[1])\n\ntotal = df.shape[0]\nfor k, p in enumerate(ax[0].patches):\n height = p.get_height()\n ax[0].text(p.get_x()+p.get_width()/2.,\n height + 3,\n \"{:1.1f}%\".format(100*height/total),\n ha=\"center\",\n fontsize=12, rotation=0)\n \nfor k, p in enumerate(ax[1].patches):\n height = p.get_height()\n ax[1].text(p.get_x()+p.get_width()/2.,\n height + 3,\n \"{:1.1f}%\".format(100*height/total),\n ha=\"center\",\n fontsize=12, rotation=70)\n \nfor label in ax[0].get_xticklabels():\n label.set_rotation(0)\nfor label in ax[1].get_xticklabels():\n label.set_rotation(0)\n\n\nax[0].set_ylabel(\"count\")\nax[0].set_title(\"ISUP Grade/Data Provider\\n\")\nax[1].set_title(\"Gleason Score/Data Provider\\n\")\nfig.suptitle(\"Relative Distribution plot (count and %)\", y=1.1)\n\nfig.tight_layout();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Inferences on the Relative Distributions\n\n#### About `ISUP Score` and `Data Provider`\n\n* In isup_grade category 0 and 1 most of the data is provided by `karolinska`.\n* In isup_grade category 3,4 and 5 most of the data is provided by `radbound`.\n\n#### About`Gleason Score` and `Data Provider`\n\n* In gleason_score category (0+0), all the data is provided by `karolinska`.\n* In gleason_score category (negative), all the data is provided by `radbound`.\n* Also in gleason_score category (3+3), karolinska is major data provider.\n* On the other hand radbound is major data provider for (4+4), (4+3), (4+5), (5+4), (5+5), (5+3), (3+5).\n\n\n# Some Advanced EDA\n\nIn this section, we are going to load some few images samples from the training image directory and learn some of their properties. The PANDA dataset are whole slide images (WSIs) stored in **Tag Image File Format** (`.tiff`). Let us recall `.tiff` is a common format for exchanging raster graphics (bitmap) images between application programs, including especially those used for scanner images. Thus, `.tiff` files can be in any of several classes, including gray scale, color palette, or RGB full color. The loading and processing of these WSIs will be performed using both `skimage` and `openslide` libraries in Python. Both libraries have some amazing and independent functions to process WSIs, for example, one of the key benefits of `openslide` is that we can load arbitrary regions of the slide, without having to load the whole image into memory. Click on [OpenSlide Python Documentation](https://openslide.org/api/python/) to read more about the OpenSlide Python. \n \n\n## Loading and Quickly Displaying of some few WSIs\n\nLet's randomly select some few WSIs from the training samples and play with them to get some few insights.","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"### Randomly sample 9 training samples from our dataframe train_ids\nnp.random.seed(13)\nWSI9 = train_ids.sample(n=9)\ndisplay(WSI9.style.background_gradient(cmap=\"Blues\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Checking of some Generalizable Properties on the Data","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"### Let's open two files from different data provider, from the 7 selected images above\nprint(f\"Data provider: {WSI9.data_provider.tolist()[0]}\\n\")\nfile_path1 = os.path.sep.join([train_dir, WSI9.index[0]+\".tiff\"]) # Full file directory\nexample_slide1 = OS.OpenSlide(file_path1) # Openining without reading the image into memory\nprint(f\"Level dimensions: {example_slide1.level_dimensions}\\n\")\n\nfor prop in example_slide1.properties.keys():\n print(\"{}: -> {}\".format(prop, example_slide1.properties[prop]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"print(f\"Data provider: {WSI9.data_provider.tolist()[1]}\\n\")\nfile_path2 = os.path.sep.join([train_dir, WSI9.index[1]+\".tiff\"]) # Full file directory\nexample_slide2 = OS.OpenSlide(file_path2) # Openining without reading the image into memory\nprint(f\"Level dimensions: {example_slide2.level_dimensions}\\n\")\n\nfor prop in example_slide2.properties.keys():\n print(\"{}: -> {}\".format(prop, example_slide2.properties[prop]))\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Some Few Insights\n\nThe outputs of the cell above show the properties of two WSIs from the 9 randomly selected training images from different data provider. What we learned from these outputs are the following:\n\n - The image dimensions are quite large, ranging typically between about 4 000 and 40 000 pixels in both x and y.\n - Each WSI has 3 levels that can be loaded, corresponding to a downsampling of about 1, 4 and 16 as shown in the output cell above. So, intermediate levels can be created by simply downsampling a higher resolution level.\n - The dimensions of each level differ based on the dimensions of the original image.\n - Biopsies can be in different rotations. This rotation has no clinical value, and is only dependent on how the biopsy was collected in the lab.\n - There are noticable color differences between the biopsies, this is very common within pathology and is caused by different laboratory procedures.\n \n\n Besside these WSI properties, I associated a $512\\times 512$ zoomed part of the image around $(x,y)=(1024,1024)$. And as you can observe on the figure below, starting from $(x,y)=(1024,1024)$ on some images, a zomm of $512\\times 512$ doesn't capture the biopsy.","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"### Visualization based on OpenSlide and Matplolib packages\ndef plot_WSIs(slides): \n fig, ax = plt.subplots(nrows=3, ncols=3, figsize=(18,18))\n for i, slide in enumerate(slides):\n file_path = os.path.sep.join([train_dir, slide+\".tiff\"]) # Full file directory\n image = OS.OpenSlide(file_path) # Openining without reading the image into memory\n \n # Creation of the patch to plot\n patch = image.read_region(location=(0,0),\n level=image.level_count-1, # Get the last level/slide\n size=image.level_dimensions[-1]) # Get the dimension corresponding of the last level\n \n # Plot the patch\n ax[i//3, i%3].imshow(patch)\n image.close()\n ax[i//3, i%3].axis(\"on\")\n \n image_id = slide\n data_provider = train_ids.loc[slide, \"data_provider\"]\n isup_grade = train_ids.loc[slide, 'isup_grade']\n gleason_score = train_ids.loc[slide, 'gleason_score']\n ax[i//3, i%3].set_title(f\"\\nID: ~{image_id[:7]}, Source: {data_provider}\\nISUP: {isup_grade}, Gleason: {gleason_score}\")\n\n fig.tight_layout()\n fig.suptitle(\"\")\n plt.show()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"plot_WSIs(WSI9.index)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Let's Resize and Visualize the 9 Randomly Selected Biopsy Images\n\n\n`skimage.io.MultiImage()` and `cv2.resize()` can be used as well instead of `openslide.OpenSlide()` and `cv2.resize()`, but the former currently fails fails on the latest version of the `Kaggle Docker image` to the following packages `imagecodecs` and `tifffile`.","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"### Visualization based on Skimage and Matplotlib packages\ndef plot_resized_biopsy(slides): \n fig, ax = plt.subplots(nrows=3, ncols=3, figsize=(18,18))\n for i, slide in enumerate(slides):\n file_path = os.path.sep.join([train_dir, slide+\".tiff\"]) # Full file directory\n image = OS.OpenSlide(file_path) # Openining without reading the image into memory\n \n # Creation of the patch to plot\n patch = image.read_region(location=(0,0),\n level=image.level_count-1, # Get the last level/slide\n size=image.level_dimensions[-1]) # Get the dimension corresponding of the last level\n \n # Resize the image patch\n image = cv2.resize(np.asarray(patch), (512, 512))\n \n # Plot the resized image patch\n ax[i//3, i%3].imshow(image) \n ax[i//3, i%3].axis(\"on\")\n \n image_id = slide\n data_provider = train_ids.loc[slide, \"data_provider\"]\n isup_grade = train_ids.loc[slide, 'isup_grade']\n gleason_score = train_ids.loc[slide, 'gleason_score']\n ax[i//3, i%3].set_title(f\"\\nID: ~{image_id[:7]}, Source: {data_provider}\\nISUP: {isup_grade}, Gleason: {gleason_score}\")\n\n fig.tight_layout()\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"plot_resized_biopsy(WSI9.index)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n### Exploration and Analysis of the Training Label Masks\n\n#### Overview\nApart from the slide-level label (present in the csv file), almost all slides in the training set have an associated mask with additional label information. The importance of these masks is they directly indicate which parts of the tissue are healthy and which are cancerous. However, the information in the masks differ from the two centers. Here is how:\n- **Radboudumc**: the prostate glands are individually labelled and the valid values are:\n - 0: background (non tissue) or unknown\n - 1: stroma (connective tissue, non-epithelium tissue)\n - 2: healthy (benign) epithelium\"\n - 3: cancerous epithelium (Gleason 3)\n - 4: cancerous epithelium (Gleason 4)\n - 5: cancerous epithelium (Gleason 5)\n- **Karolinska**: the regions are labelled and the valid values are:\n - 0: background (non tissue) or unknown\n - 1: benign tissue (stroma and epithelium combined)\n - 2: cancerous tissue (stroma and epithelium combined)\n\nThe label masks provided by Radboudumc center were semi-automatically generated by several deep learning algorithms, contain noise, and can be considered as weakly-supervised labels. The label masks provided by Karolinska center were semi-autotomatically generated based on annotations by a pathologist.\n\nThe label masks are stored in an RGB format so that they can be easily opened by image readers. The label information is stored in the red (R) channel, the other channels are set to zero and can be ignored. As with the slides itself, the label masks can be opened using OpenSlide.\n\n#### Loading of label masks and Visualization","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"def plot_biopsy_masks(slides): \n fig, ax = plt.subplots(3,3, figsize=(18,18))\n for i, slide in enumerate(slides):\n \n file_path = os.path.sep.join([train_mask_dir, slide+\"_mask.tiff\"]) # Full file directory\n biopsy_mask = OS.OpenSlide(file_path) # Openining without reading the image into memory\n \n # Creation of the patch to plot\n mask_data = biopsy_mask.read_region(location=(0,0),\n level=biopsy_mask.level_count - 1, # Get the last level/slide\n size=biopsy_mask.level_dimensions[-1]) # Get the dimension corresponding of the last level\n \n # Plot\n cmap = mpl.colors.ListedColormap([\"black\", \"gray\", \"green\", \"yellow\", \"orange\", \"red\"])\n ax[i//3, i%3].imshow(np.asarray(mask_data)[:,:,0], cmap=cmap, interpolation=\"nearest\", vmin=0, vmax=5) \n biopsy_mask.close() \n ax[i//3, i%3].axis(\"on\")\n \n image_id = slide\n data_provider = train_ids.loc[slide, \"data_provider\"]\n isup_grade = train_ids.loc[slide, \"isup_grade\"]\n gleason_score = train_ids.loc[slide, \"gleason_score\"]\n ax[i//3, i%3].set_title(f\"\\nID: {image_id[:7]} Source: {data_provider}\\nISUP: {isup_grade} Gleason: {gleason_score}\")\n \n fig.tight_layout() \n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"plot_biopsy_masks(WSI9.index)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Prepare the Training Sample Images with their Corresponding Masks.\n\nRecall we showed early that there are exactly 100 training sample images whose the training label masks are not available in the mask directory. So, we are going rearrange each training sample with its corresponding mask.","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"train_df = train_ids.reset_index() # In train_ids, I set index to \"image_id\"\n\nmasks = os.listdir(train_mask_dir)\nmasks_df = pd.DataFrame(data={\"mask_id\": masks})\nmasks_df[\"image_id\"] = masks_df.mask_id.apply(lambda x: x.split(\"_\")[0]) # Recall mask_id=image_id+\"_mask.tiff\"\n\ntrain_df = pd.merge(train_df, masks_df, on=\"image_id\", how=\"outer\")\nprint(f\"We have {train_df.shape[0]} training sample images and {masks_df.shape[0]} masks. So, there will be exactly {len(train_df[~train_df.mask_id.isna()])} images in the final training samples.\")\ndisplay(train_df.head(10).style.background_gradient(cmap=\"Blues\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Plotting of some Image and Masks According to `sup_grade` Category\n\n\n**Auxiliary functions for this task:**","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"def load_and_resize_biopsy(img_id):\n \n file_path = os.path.sep.join([train_dir, img_id+\".tiff\"]) # Full file directory\n biopsy_img = OS.OpenSlide(file_path) # Openining without reading the image into memory\n\n # Creation of the patch to plot\n patch = biopsy_img.read_region(location=(0,0),\n level=biopsy_img.level_count-1, # Get the last level/slide\n size=biopsy_img.level_dimensions[-1]) # Get the dimension corresponding of the last level\n\n # Resize the image patch\n image = cv2.resize(np.asarray(patch), (512, 512))\n \n return image\n\ndef load_and_resize_biopsy_mask(img_id):\n \n file_path = os.path.sep.join([train_mask_dir, img_id+\"_mask.tiff\"]) # Full file directory\n biopsy_mask = OS.OpenSlide(file_path) # Openining without reading the image into memory\n\n # Creation of the patch to plot\n patch = biopsy_mask.read_region(location=(0,0),\n level=biopsy_mask.level_count-1, # Get the last level/slide\n size=biopsy_mask.level_dimensions[-1]) # Get the dimension corresponding of the last level\n\n # Resize the mask patch\n mask = cv2.resize(np.asarray(patch), (512, 512))[:,:,0]\n \n return mask","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Visualization","execution_count":null},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"### Visualization\ndata_providers = train_df.data_provider.unique().tolist()\n# cmap = mpl.colors.ListedColormap([\"black\", \"gray\", \"green\", \"yellow\", \"orange\", \"red\"])\ncmap_rad = mpl.colors.ListedColormap([\"white\", \"lightgrey\", \"green\", \"orange\", \"red\", \"darkred\"])\ncmap_kar = mpl.colors.ListedColormap([\"white\", \"green\", \"red\"])\nlabels = []\nfor grade in range(train_ids.isup_grade.nunique()):\n fig, ax = plt.subplots(nrows=4, ncols=4, figsize=(22, 22))\n\n for i, row in enumerate(ax):\n idx = i//2\n temp_idx = (train_df.isup_grade == grade) & (train_df.data_provider == data_providers[idx])\n temp = train_df[temp_idx].image_id.tail(4).reset_index(drop=True)\n if i%2 < 1:\n labels.append(f\"{data_providers[idx]}\\n(image)\")\n for j, col in enumerate(row):\n col.imshow(load_and_resize_biopsy(temp[j]))\n col.set_title(f\"\\nID: {temp[j][:13]} $\\cdots$\")\n \n else:\n labels.append(f\"{data_providers[idx]}\\n(mask)\")\n for j, col in enumerate(row):\n if data_providers[idx] == \"radboud\":\n col.imshow(load_and_resize_biopsy_mask(temp[j]), \n cmap = cmap_rad, \n norm = mpl.colors.Normalize(vmin=0, vmax=5, clip=True))\n else:\n col.imshow(load_and_resize_biopsy_mask(temp[j]),\n cmap = cmap_kar,\n norm = mpl.colors.Normalize(vmin=0, vmax=2, clip=True))\n \n gleason_score = train_ids.loc[temp[j], \"gleason_score\"]\n col.set_title(f\"\\nID: {temp[j][:13]} $\\cdots$\")\n \n for row, r in zip(ax[:,0], labels):\n row.set_ylabel(r, rotation=0, size=\"large\", labelpad=30, fontsize=20)\n \n fig.tight_layout()\n fig.suptitle(f\"ISUP Grade {grade}\", y=1.01, fontsize=23)\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Overlaying of Training Label Masks on their Corresponding Whole Slide Images\n\nGiven the knowledge of the importance of masks and also the fact that they have the same dimension as their corresponding slides we can overlay them on the tissue to directly see areas which are cancerous. This overlay can help you identifying the different growth patterns. To do this, we load both the mask and the biopsy with `OS.OnpenSlide()` and merge them using `PIL.Image()`.\n\n#### Auxiliary functions for this tasks\n\nAfter creating the mask patch to visualize, we will apply `mask.split()` function that will split the image, `mask`, into individual bands/channels. This method returns a tuple of individual image bands from an image. For example, splitting an \"RGB\" image creates three new images each containing a copy of one of the original bands/channels (red, green, blue).","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"def mask_on_slide_overlayer(images, center=\"radboud\", fig_title=\"\", alpha=0.8, max_size=(800, 800)):\n \"\"\"Show a mask overlayed on a slide.\"\"\"\n fig, ax = plt.subplots(nrows=3,ncols=3, figsize=(18,18))\n \n \n for i, image_id in enumerate(images):\n # Open a slide/wsi\n wsi_file_path = os.path.sep.join([train_dir, image_id+\".tiff\"]) # Full file directory\n biopsy_img = OS.OpenSlide(wsi_file_path) # Openining without reading the image into memory\n # Open the corresponding mask\n mask_file_path = os.path.sep.join([train_mask_dir, image_id+\"_mask.tiff\"])\n biopsy_mask = OS.OpenSlide(mask_file_path)\n\n # Creation of the patch to visualize\n patch_img = biopsy_img.read_region(location=(0,0),\n level=biopsy_img.level_count-1, # Get the last level/slide\n size=biopsy_img.level_dimensions[-1]) # Get the dimension corresponding of the last level\n \n patch_mask = biopsy_mask.read_region(location=(0,0),\n level=biopsy_mask.level_count-1,\n size=biopsy_mask.level_dimensions[-1])\n\n\n # Split the patch mask into channels\n patch_mask = patch_mask.split()[0]\n \n # Create alpha mask\n alpha_int = int(round(255*alpha))\n if center == \"radboud\":\n alpha_content = np.less(patch_mask.split()[0], 2).astype('uint8') * alpha_int + (255 - alpha_int)\n elif center == \"karolinska\":\n alpha_content = np.less(patch_mask.split()[0], 1).astype('uint8') * alpha_int + (255 - alpha_int)\n\n alpha_content = Image.fromarray(alpha_content)\n preview_palette = np.zeros(shape=768, dtype=int)\n\n if center == \"radboud\":\n # Mapping: {0: background, 1: stroma, 2: benign epithelium, 3: Gleason 3, 4: Gleason 4, 5: Gleason 5}\n preview_palette[0:18] = (np.array([0, 0, 0, 0.5, 0.5, 0.5, 0, 1, 0, 1, 1, 0.7, 1, 0.5, 0, 1, 0, 0]) * 255).astype(int)\n elif center == \"karolinska\":\n # Mapping: {0: background, 1: benign, 2: cancer}\n preview_palette[0:9] = (np.array([0, 0, 0, 0, 1, 0, 1, 0, 0]) * 255).astype(int)\n\n patch_mask.putpalette(data=preview_palette.tolist())\n mask_rgb = patch_mask.convert(mode=\"RGB\")\n \n # Overlay the mask on its corresponding slide\n overlayed_image = Image.composite(image1=patch_img, image2=mask_rgb, mask=alpha_content)\n overlayed_image.thumbnail(size=max_size, resample=0)\n\n # Plot the overlayed image\n ax[i//3, i%3].imshow(overlayed_image) \n biopsy_img.close()\n biopsy_mask.close() \n ax[i//3, i%3].axis(\"on\")\n \n data_provider = train_ids.loc[image_id, \"data_provider\"]\n isup_grade = train_ids.loc[image_id, \"isup_grade\"]\n gleason_score = train_ids.loc[image_id, \"gleason_score\"]\n ax[i//3, i%3].set_title(f\"\\nID: {image_id[:7]} $\\cdots$, Source: {data_provider}\\nISUP: {isup_grade} Gleason: {gleason_score}\")\n \n fig.suptitle(fig_title, y=1.01, fontsize=23)\n fig.tight_layout()\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Visualization","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"fig_title = \"Some Examples of Overlayed Images\"\nmask_on_slide_overlayer(images=WSI9.index, fig_title=fig_title)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Some Few Insights\n\nAs you can observe on the outputs of the cell above, overlaying mask on slide, there are some few pen markings on the slide (dark green smudges). As explained in the overview on competition, these markings are not part of the tissue but were made by the pathologists who originally checked this case. As the competition organizers described, *slightly different procedures were in place for the images used in the test set than the training set. Some of the training set images have stray pen marks on them, but the test set slides are free of pen marks., these pen markings are available on some slides in the training set.*\n\n### Let's Explore few of these Images with Pen Markers\n","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"WSI9.index[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"pen_marked_images = [\"ca0798453868081bc8aeeabb01847d4e\",\n \"ff10f937c3d52eff6ad4dd733f2bc3ac\",\n \"e9a4f528b33479412ee019e155e1a197\",\n \"fd6fe1a3985b17d067f2cb4d5bc1e6e1\",\n \"f39bf22d9a2f313425ee201932bac91a\",\n \"fb01a0a69517bb47d7f4699b6217f69d\",\n \"ebb6a080d72e09f6481721ef9f88c472\",\n \"feee2e895355a921f2b75b54debad328\",\n \"ebb6d5ca45942536f78beb451ee43cc4\"]\n\n\nfig_title = \"Some Examples of Pen Marked Images\"\nmask_on_slide_overlayer(images=pen_marked_images, fig_title=fig_title)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Distributions of the Dimensions (`width`,`hight`) of Training Samples\n\n### Extraction of Dimensions (`width`,`hight`) and Merging with the Train Data","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"trn_df = train_df.copy()\ndims, spacings = [], []\n\nfor img_id in trn_df.reset_index().image_id:\n # Open a slide/wsi\n wsi_file_path = os.path.sep.join([train_dir, img_id+\".tiff\"]) # Full file directory\n biopsy_img = OS.OpenSlide(wsi_file_path) # Openining without reading the image into memory\n \n spacing = 1 / (float(biopsy_img.properties[\"tiff.XResolution\"]) / 10000)\n dims.append(biopsy_img.dimensions)\n spacings.append(spacing)\n biopsy_img.close()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"trn_df[\"spacing\"] = spacings\ntrn_df[\"width\"] = [i[0] for i in dims]\ntrn_df[\"height\"] = [i[1] for i in dims]\n\ndisplay(trn_df.head(10).style.background_gradient(cmap=\"Blues\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Visualization of some Insighful Distributions\n\n\n#### Auxiliary Function","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"def plot_distribution_grouped(feature, feature_group, ax):\n for feat in trn_df[feature_group].unique():\n df = trn_df.loc[trn_df[feature_group] == feat]\n sns.kdeplot(df[feature], label=feat, ax=ax, shade=True)\n ax.set_title(f\"Images {feature}\\ngrouped by {feature_group}\\n\")\n ax.legend()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(15,6), dpi=200, sharey=True)\n###\nsns.kdeplot(trn_df[\"width\"], ax=ax[0], shade=True, label=\"width\")\nsns.kdeplot(trn_df[\"height\"], ax=ax[0], shade=True, label=\"height\")\nax[0].set_xlabel(\"dimension\")\nax[0].set_title(\"Images Width and Height\\n\")\nax[0].legend()\n###\nplot_distribution_grouped(feature=\"width\", feature_group=\"data_provider\", ax=ax[1])\nplot_distribution_grouped(feature=\"height\", feature_group=\"data_provider\", ax=ax[2])\n\nfig.suptitle(\"Distribution Plots\", y=1.1)\nfig.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(12,6), dpi=200, sharey=True)\n###\nplot_distribution_grouped(feature=\"width\", feature_group=\"isup_grade\", ax=ax[0])\nplot_distribution_grouped(feature=\"height\", feature_group=\"isup_grade\", ax=ax[1])\n\nfig.suptitle(\"Distribution by ISUP Grade\", y=1.1)\nfig.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(12,6), dpi=200, sharey=True)\n###\nplot_distribution_grouped(feature=\"width\", feature_group=\"gleason_score\", ax=ax[0])\nplot_distribution_grouped(feature=\"height\", feature_group=\"gleason_score\", ax=ax[1])\n\nfig.suptitle(\"Distribution by Gleason Score\", y=1.1)\nfig.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Inferences on the Distributions\n\nWe can see that overall, as well as at `ISUP Grade` and `Gleason Socre` levels, width and height have almost similar distribution.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Inspirational References\n\n\nThe following kernels were of great help in the creation of this notebook:\n\n1. The main domain knowledge info has been taken from the following sources:\n - https://www.kaggle.com/c/prostate-cancer-grade-assessment/overview\n - https://www.kaggle.com/c/prostate-cancer-grade-assessment/data\n - https://zenodo.org/record/3715938#.XxR-85MzZQL\n - https://emedicine.medscape.com/article/1612022-overview\n1. https://www.kaggle.com/gpreda/panda-challenge-starting-eda\n1. https://www.kaggle.com/iamleonie/panda-eda-visualizations-suspicious-data\n ","execution_count":null},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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":4,"nbformat_minor":4}