phenotypes_eda.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "997e38d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "library(data.table)\n",
    "library(xml2)\n",
    "library(XML)\n",
    "library(dplyr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad5d9357",
   "metadata": {},
   "source": [
    "# Read GTEx V8 Subject Phenotype Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd4d4c99",
   "metadata": {},
   "outputs": [],
   "source": [
    "phenos <- fread(\"data/phenotypes/GTEx_Subject_Phenotypes.GRU.txt.gz\"\n",
    "                ,skip=\"dbGaP_Subject_ID\"\n",
    "               ,header=TRUE)\n",
    "head(phenos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19027819",
   "metadata": {},
   "outputs": [],
   "source": [
    "with(phenos, table(SEX))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a1657d1",
   "metadata": {},
   "source": [
    "Males are marked as \"1\" and females as \"2\" (according to dictionary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a247d4bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "males <- subset(phenos, phenos$SEX %in% 1)\n",
    "females <- subset(phenos, phenos$SEX %in% 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a0f1464",
   "metadata": {},
   "source": [
    "# Load Sample Attributes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3d322e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "sampleattr <- fread(\"data/phenotypes//GTEx_Sample_Attributes.GRU.txt.gz\"\n",
    "                   ,skip=\"dbGaP_Sample_ID\"\n",
    "                   ,header=TRUE)\n",
    "head(sampleattr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "664f51b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "sampleattrxml <- read_xml(\"data/phenotypes//GTEx_Sample_Attributes.data_dict.xml\")\n",
    "sampleattrxml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "697922bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "samplecolnames <- xml_text(xml_find_all(sampleattrxml, xpath=\"//name\"))\n",
    "sampledesc <- xml_text(xml_find_all(sampleattrxml, xpath=\"//description\")[-1])\n",
    "sampleattr_tb <- tibble(attrname = samplecolnames, description=sampledesc)\n",
    "sampleattr_tb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d70bc03",
   "metadata": {},
   "outputs": [],
   "source": [
    "sampleattr[,c('dbGaP_Sample_ID','SAMPID','SMSMPSTE','SMTS','SMTSD')]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0008fd6b",
   "metadata": {},
   "source": [
    "So the variable for the tissue type seems to be `SMTSD`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84e81bbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "table(sampleattr$SMTSD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3accfd80",
   "metadata": {},
   "outputs": [],
   "source": [
    "brainsamples <- subset(sampleattr, sampleattr$SMTS %in% c(\"Brain\",\"Nerve\"))\n",
    "table(brainsamples$SMTSD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7af5429e",
   "metadata": {},
   "outputs": [],
   "source": [
    "head(brainsamples[,c('dbGaP_Sample_ID','SAMPID','SMSMPSTE','SMTS','SMTSD')])\n",
    "paste(\"Total number of brain and nerve samples:\", nrow(brainsamples))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49428cd9",
   "metadata": {},
   "source": [
    "# Obtain the Subject and Sample IDs Needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48613ef3",
   "metadata": {},
   "outputs": [],
   "source": [
    "subjid <- gsub(\"(GTEX-.*)-(.*)*\",\"\\\\1\",brainsamples$SAMPID)\n",
    "paste(\"Number of subjects: \", length(unique(subjid)))\n",
    "x <- merge(tibble(SUBJID=subjid), phenos, by=\"SUBJID\")\n",
    "#head(x[,c(\"SUBJID\",\"SEX\",\"AGE\",\"RACE\",\"ETHNCTY\")])\n",
    "x2 <- bind_cols(brainsamples, x)\n",
    "\n",
    "tissue <- x2$SMTSD\n",
    "tissue <- gsub(\"-\",\"\",tissue)\n",
    "tissue <- gsub(\"[[:space:]]\", \".\", tissue)\n",
    "tissue <- gsub(\"\\\\.\\\\.?\",\".\",tissue)\n",
    "tissue <- gsub(\"\\\\(|\\\\)\",\"\",tissue)\n",
    "sex <- ifelse(x2$SEX==1,\"male\",\"female\")\n",
    "outdir=\"data/phenotypes/sample_subsets/\"\n",
    "fname <- paste0(outdir,tissue,\".\",sex)\n",
    "\n",
    "x3 <- bind_cols(x2, fname=fname)\n",
    "\n",
    "colist <- c('dbGaP_Sample_ID','SAMPID','SUBJID',\"SEX\",\"AGE\",\"RACE\",\"ETHNCTY\",'SMSMPSTE','SMTS','SMTSD','fname')\n",
    "\n",
    "head(x3[, ..colist])\n",
    "nrow(x3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d694eb7",
   "metadata": {},
   "source": [
    "Have to get separate sample lists for each brain/nerve tissue types and also separate out by gender.  \n",
    "In addition, need to ensure sample lists have both genotypes and expression data.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e841b3ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "famfile <- \"data/genotypes/phg001219.v1.GTEx_v8_WGS.genotype-calls-vcf.c1/GTEx_Analysis_2017-06-05_v8_WholeGenomeSeq_838Indiv_Analysis_Freeze.SHAPEIT2_phased_bplink.fam\"\n",
    "fam <- fread(famfile)\n",
    "head(fam)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44b6c176",
   "metadata": {},
   "outputs": [],
   "source": [
    "#by_tissue_sex <- x3 %>% group_by(SMTSD,SEX)\n",
    "#bytissue %>% summarise(n=n())\n",
    "fnames <- unique(x3$fname)\n",
    "for(i in fnames) {\n",
    "    s <- subset(x3, x3$fname %in% i)\n",
    "    s <- merge(s,fam[,2],by.x=\"SUBJID\",by.y=\"V2\")\n",
    "    if(nrow(s)>0) {\n",
    "        fwrite(s[,'SAMPID'], paste0(i,\".SAMPID_list.txt\")\n",
    "               ,quote=F, row.names=F, col.names=F, sep=\"\\t\")\n",
    "        fwrite(s[,c('SAMPID','SUBJID')], paste0(i,\".SAMPID_lookup.txt\")\n",
    "               ,quote=F, row.names=F, col.names=T, sep=\"\\t\")\n",
    "        fwrite(s[, ..colist], paste0(i,\".sample_attributes.txt\")\n",
    "              ,quote=F, row.names=F, col.names=T, sep=\"\\t\")\n",
    "    }\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11e5ab88",
   "metadata": {},
   "source": [
    "# Subject Phenotype Data Dictionary"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21adb049",
   "metadata": {},
   "source": [
    "It is best to go to the following link:  \n",
    "https://ftp.ncbi.nlm.nih.gov/dbgap/studies/phs000424/phs000424.v8.p2/pheno_variable_summaries/phs000424.v8.pht002742.v8.GTEx_Subject_Phenotypes.data_dict.xml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66cf5808",
   "metadata": {},
   "outputs": [],
   "source": [
    "phenoxml <- read_xml(\"data/phenotypes//GTEx_Subject_Phenotypes.var_report.xml\")\n",
    "phenoxml_variable <- xml_find_all(phenoxml, xpath=\"//variable\")\n",
    "\n",
    "#phenoxml_variable_text <- xml_text(phenoxml_variable)\n",
    "#phenoxml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1e6a6a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#phenolist <- read_xml(\"data/phenotypes//GTEx_Subject_Phenotypes.var_report.xml\") %>% as_list()\n",
    "#phenolist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "668e5a3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#phenoxml_variables <- xmlToDataFrame(nodes=getNodeSet(phenoxml, \"//variable\"))\n",
    "#phenoxml_variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "620fcff9",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}