Lesson goal: take your dataset and add columns from another datasets. You will add columns by matching shared values between the two datasets (usually columns with gene names or Ensembl Gene IDs).
Download the following files to folder "data":
- "Additional File 2: Table S1" Excel file from Gonzalez et al. "Sex differences in the late first trimester human placenta transcriptome", Biology of Sex Differences (2018). It is under section "Additional files".
- "proteinatlas.tsv.zip" from The Human Protein Atlas downloadable data page.
- Save the Human Protein Atlas page as a pdf (e.g. print to pdf from Chrome browser) or manually record the database version that you are using. It is good practice to document the database version you are using for later citation and methods sections.
Save the files to a data folder that you can access through R Studio.
In my examples below, my working directory will have subfolders named "data" (for my data) and "hpa2026-07" (for the month when I downloaded Human Protein Atlas files). You can do the same, or leave everything in your working directory without using subfolders.
Set your working directory (the folder where your want R to look for your data):
setwd("C:\Dropbox\POST-DOC\Bioinfo\HPA") ## backward slashes will give you an errorsetwd("C:/Dropbox/POST-DOC/Bioinfo/HPA") ## forward slashes are correctClean up your starting dataset for R input:
To feed data tables into R, the tables need to be perfectly rectangular. No extra rows or columns below or to the left/right of the table are allowed.
Additional File 2: Table S1 has extra rows at the bottom that aren't part of the table. Please select, right-click, and delete these rows. Extra rows above are ok and can be skipped during data input. For additional details and screenshots, see lesson #1 section "Clean Additional File 2 for use with R programming".
 |
| Extra rows below the table (highlighted red for illustrative purposes) will be deleted. They will interfere with R input. Delete the whole rows, not just the values they contain. |
Lesson #1 taught you to save as a comma-separated (.csv) file, a very common and simple file type for bioinformatics data. If you already loaded this table as variable df from lesson #1, you can skip this section use that df.
Otherwise, save the updated table as Excel file "gonzalez2018-tableS1.xlsx" and load Excel data using openxlsx R package instead. This code below will only work to open Excel files with file extension .xlsx but not .csv or .xls or .xlm or other file extensions.
## Install and load package to read Excel files -----------------install.packages("openxlsx") # just do oncelibrary(openxlsx)
## Input my data ------------------------------------------------
mytablefile <- "./data/gonzalez2018-tableS1.xlsx" ## example path to my data file, but use your own xlsx path. Remember that "./" means start at your current working directory.
file.exists(mytablefile) ## If TRUE, R can find your file (good, continue).
## If FALSE, look for typos or path mistakes.df <- read.xlsx(xlsxFile = mytablefile, sheet=1, startRow=3, na.strings = c("NA", "", " "))colnames(df)
If your file.exists(mytablefile) results in FALSE, rule out these mistakes:
- Typo in the file name
- Typo in the file extension (are you using the previous .csv file instead of an .xlsx file?)
- Working directory or subfolder not available. Check what R can see using: dir() to see the contents of your working directory. Check dir("data") to see the contents of subfolder "data", or whatever your subfolder might be named.
- Working directory not what you expected. Check your working directory: getwd()
Variable sheet=1 tells function open.xlsx which Excel sheet to load, in case there are multiple sheets.
Variable startRow=3 tells function open.xlsx which row contains our column titles and where our table begins. This is similar to lesson #1's variable skip in function read.csv. Each function needs to know where your table begins, but they define the start slightly differently.
 |
| Red arrow shows that the header row (column titles) are on row 3. |
Variable na.strings tells R how the Excel file denotes missing values. In this example, any cell that is just the value "NA", nothing (""), or a space (" ") will be called a missing value.
After successfully loading data into variable df, check that your column names are as expected with:
colnames(df)
If your column names are not what you expected, edit the value of startRow and try again.
Load packages and unzip the Human Protein Atlas data
Because the Human Protein Atlas supplies files as zip-compressed tables, you need to first unzip each file. You can do this manually with free software such as 7Zip or PeaZip, or with R code as shown below:
## Install packages (just do once) -------------------------------## remove first hashtag and install if needed
#install.packages("tidyverse")
#install.packages("data.table")
#install.packages("HGNChelper")
## Load packages (do every time you run this code) ---------------
library(openxlsx)library(tidyverse)
library(data.table)
library(HGNChelper)
#### SET WORKING DIRECTORY ####################################################
## be in directory where folder where zip files are located
setwd("../HPA/hpa2026-07/") ## EDIT THIS LINE FOR YOUR OWN PATH WITH HPA FILES
dir() ## see what files are available
#### UNZIP FILES FIRST TIME, ALL AT ONCE ######################################
## make temporary directory for unzipped file contents
tempdir <- "./temp/"
dir.create(tempdir)
## identify zip files by looking for ".zip" at the end of each directory item
allzip <- unlist(grep("[.]zip$", dir(), value=TRUE))
allzip
## unzip all files identified above
for (z in allzip) {
print(z)
print(unzip(z, exdir = tempdir, overwrite=FALSE, list=TRUE)) ## doesn't unzip, only lists files
unzip(z, exdir = tempdir, overwrite=TRUE) ## actually extract contents
cat("\n-----------\n")
}
#### LOAD DATA FROM FILES #####################################################
## Human protein atlas data
hpa <- read.table("./temp/proteinatlas.tsv",
sep="\t", ## tab-separated files
quote = "\"", ## if this isn't added, there's a row with a single ' that causes truncation
header=TRUE,
stringsAsFactors=FALSE)
colnames(hpa)[1:30] ## see the first 30 column names
dim(hpa) ## how many rows and columns?
## Protein expression (IHC)
ihc <- read.table("./temp/normal_ihc_data.tsv",
sep="\t", ## tab-separated files quote = "\"", ## if this isn't added, there's a row with a single ' that causes truncation
header=TRUE,
stringsAsFactors=FALSE)
head(ihc, 10) ## view the first 10 rows
table(ihc$Level, ihc$Tissue) ## make a table to count values from Level and Tissue columns
Summarize protein classes and secretome locations from Human Protein Atlas
This HPA dataset contains a column for all protein classes that correspond to each protein-coding gene (each row), but the protein classes are listed together as comma-separated values. That makes it hard to view and use. Below, make csv spreadsheets to separate each row's values by the comma delimiter (sep=",") and count each time a value appears in the dataset.
Do the same thing for the secretome location column.
Code below:
#### EXAMINE DATA CATEGORIES ##################################################
## Count protein classes ---------------------------------
head(hpa$Protein.class, 20)
protein_class_counts <- hpa %>%
select(Protein.class) %>%
filter(!is.na(Protein.class), Protein.class != "") %>%
separate_rows(Protein.class, sep = ",") %>%
mutate(Protein.class = str_trim(Protein.class)) %>%
count(Protein.class, sort = TRUE) %>%
mutate(percent = scales::percent(n / sum(n), accuracy = 0.1))
## output as spreadsheet
write.csv(protein_class_counts, "summary_Protein.class_counts.csv")
## Examine secretome location ----------------------------
head(hpa$Secretome.location, 20)
secretome_location_counts <- hpa %>%
select(Secretome.location) %>%
filter(!is.na(Secretome.location), Secretome.location != "") %>%
separate_rows(Secretome.location, sep = ",") %>%
mutate(Secretome.location = str_trim(Secretome.location)) %>%
count(Secretome.location, sort = TRUE) %>%
mutate(percent = scales::percent(n / sum(n), accuracy = 0.1))
## output as spreadsheet
write.csv(secretome_location_counts, "summary_Secretome.location_counts.csv")
New tidyverse syntax! Here you are using a new R syntax called "piping" which uses %>% to feed the result of each function into the next function. This syntax is not part of base R programming (base meaning default), but it is introduced by the R tidyverse collection of packages. These are common packages for R programming that supplement the base R language with additional useful tools.
If you get this error, it is because you don't have tidyverse R packages loaded in your current R session:
Solution to this error: install.packages("tidyverse") ## install now if you skipped earlier
library(tidyverse) ## load now if you didn't load it earlier
Read more about R's tidyverse:
- R for Graduate Students: Introduction to Tidyverse
- Tidyverse Skills for Data Science: Chapter 1 Introduction to the Tidyverse
- Rafa lab's Introduction to Data Science, section 4, "The tidyverse"
Create columns for specific protein classes and protein localizations
Let's say that you are interested identifying all proteins from specific protein classes, for example all "Plasma membrane" proteins. You want to be able to identify specific text strings inside the column "Protein.classes". That is what we will do below.
## PREPARE DATAFRAME FROM HPA FOR MATCHING ####################################
## Select columns to add to to your dataset. MUST also include a column to match to your dataset:
hpa.keepcol = c(
"Gene.synonym",
"Gene",
"Ensembl",
"Uniprot",
"Gene.description",
"Subcellular.location",
"Secretome.location",
"Protein.class"
)
## Check for typos: do all your column names exist in hpa? Want zero result, meaning that all columns in hpa.keepcol are in colnames(hpa):
setdiff(hpa.keepcol, colnames(hpa))
## Create a subset of hpa:
hpa_subset <- hpa[, hpa.keepcol]
## Optional: view the first 200 rows of your subset:
#View(hpa_subset[1:200,])
## Custom function ----------------------------------------------------
## Custom function to find a text patterns inside columns
## "Protein.class", "Subcellular.location", and "Secretome.location",
## then return a text pattern as a value:
grepmirror <- function(pattern, returnme) {
x <- ifelse(
## conditional statement trying to match text pattern to any
## of the following, where the bar | means "OR":
grepl(pattern, hpa_subset$Protein.class, ignore.case = TRUE) |
grepl(pattern, hpa_subset$Subcellular.location, ignore.case = TRUE) |
grepl(pattern, hpa_subset$Secretome.location, ignore.case = TRUE),
## if match found, return text string saved to variable 'returnme':
returnme,
## if no match, return an empty string:
"")
## function returns value of x
return(x)
}
## Add columns using function grepmirror(): --------------------------
hpa_subset$Membrane <- grepmirror("Predicted membrane proteins", "Membrane")
hpa_subset$Plasma_membrane <- grepmirror("Plasma membrane", "PM")
hpa_subset$Secreted <- grepmirror("secreted", "Secreted")
hpa_subset$PM_or_Secreted <- ifelse(hpa_subset$Plasma_membrane=="PM" |
hpa_subset$Secreted=="Secreted", TRUE, "")
## Optional: view the first 200 rows of your subset:
#View(hpa_subset[1:200,])
Match example #1: Select shared columns for data matching and add HPA columns to example Table S1
Continuing with Table S1, which conveniently contains Ensembl Gene IDs from RNA-seq data, we can directly match to Ensembl Gene IDs from the Human Protein Atlas dataset.
Check your column name spelling and column values:
 |
| Check the columns that you will use for matching the two datasets. Values should be similar and match-able. |
Gene symbols are acceptable for matching but beware that official gene names change, so you may have issues if the two datasets are based on very different human genome versions (e.g. 5+ years difference).
Ensembl Gene IDs (from ensembl.org) or other accession numbers are preferred because the meaning of these is more stable. By the way, the spelling is supposed to be "Ensembl" and not "Ensemble", but the Table S1 column name has a typo.
The function left_join() is a tidyverse function that adds columns from y (right) to x (left) dataframes, so which dataframes are assigned to x and y matters. The by variable indicates the column names to use for the matching. If dataframes have the same column name, you can assign it as a string variable, for example:
by="Ensembl"
But otherwise use pattern:
by = c("column_in_x" = "column_in_y")
...replacing the text inside the quotes with the specific column names in x and y (order matters).
Merge code:
## JOIN STEP! ADD SELECTED HPA COLUMNS TO YOUR DATA ###########################
## left_join is a tidyverse function that adds columns from y (right) to x (left) dataframes
## so which dataframes are assigned to x and y matters:
df_HPA <- left_join(x = df,
y = hpa_subset, by = c("Ensemble.Gene.ID" = "Ensembl"))
## Check the row numbers to make sure they look ok
## Expect no change in number of rows if your matches were 1:1
## But you might have many extra rows if you had multiple matches
dim(df)
dim(df_HPA)
dim(hpa_subset)
## Compare the number of unique Ensembl Gene IDs:
length(unique(df$Ensemble.Gene.ID))length(unique(df_HPA$Ensemble.Gene.ID))length(unique(hpa_subset$Ensembl))
## Also check your new dataset to look for any problems (e.g. empty columns indicating a matching issue):
View(df_HPA[1:100, ]) ## view first 100 rows
#View(df_HPA) # optionally view everything (commented out with a hashtag)
Review your merge results:
You started with the left dataframe (x = df), added columns from the right dataframe (y = hpa_subset), and created the new merged dataframe df_HPA by matching Ensembl Gene ID values.
Since you are using function left_join(), you expect df_HPA to be the left dataframe with a few additional columns. Your row numbers should not decrease from the left input, only stay the same or increase. The right input has many more rows but we are not using all of them so that's fine. We are only using rows that match by Ensembl Gene ID.
|
| Check dimensions (# rows, # columns) as well as unique values used for the matching. |
Also review the following:
- colnames(df_HPA) # Is the data messier now due to duplicate columns?
- View(df_HPA) # Actually look at the data to review the match results. Do some rows have only NA values (missing values) from the right dataframe? That indicates that they didn't match.
Some rows likely won't match for a good reason. For example, the left dataframe is total RNA-seq data which includes non-coding genes, and the right dataframe is from a protein database so it focuses on only coding genes. It makes sense that some non-coding genes may be missing in the Human Protein Atlas.
Alternative match by multiple columns:
Suppose that you have multiple columns that you can use for matching. If you want to match by multiple columns (requiring all values to match), this is how you do it:
df_HPA2 <- left_join(x = df,
y = hpa_subset,
by = c("Ensemble.Gene.ID" = "Ensembl",
"Gene.Symbol" = "Gene"))
The gene symbol column is named "Gene.Symbol" on my left dataframe and "Gene" on my right dataframe, so I need to tell R to treat these columns as equivalent for matching purposes. If the columns had the same name, I could just write by = "Gene" or similar.
If you want to match using only the gene symbols, ignoring Ensembl Gene IDs entirely, this is how you do it:
df_HPA3 <- left_join(x = df,
y = hpa_subset,
by = c("Gene.Symbol" = "Gene"))
Always check your merge results! Check the number of rows in the inputs and output, the number of unique values (in this case also check unique values of column "Gene.Symbol"), and the actual match results using View(df_HPA2) or View(df_HPA3) to see if you have more or less NA values.
You can also combine the table() and is.na() functions to count the number of NA values in specific columns. Here, we are adding column "Uniprot" from the right dataframe hpa_subset so checking NA values in that column for the merge result can help us identify which merge strategy produced the most matches.
I tried merging three different ways:
- Match by the Ensembl Gene ID only
- Match Ensembl Gene ID and Gene Symbol
- Match Gene Symbol only
Here are the counts of missing values for column "Uniprot". The table value TRUE indicates how many Uniprot rows are NA (missing) values in the final merge. The sum added by addmargin() shows how many total rows we have. Which merge result do you think looks better?
 |
| Comparing the results of three merge strategies with different by variable inputs. |
Less missing values is better in this example, so merging by only the Ensembl Gene ID (df_HPA) worked best. It resulted in only 3327 missing values in column "Uniprot", less than the other strategies.
Merging by only gene symbol (df_HPA3) resulted in 6 extra rows, likely because some gene symbols are not unique so the function matched multiple times.
## OUTPUT ####################################################################
## create output filename of your choice:
fn.out = "myoutput.xlsx"
## alternatively, create output filename from your input filename,
## but replace ".xlsx" with "_ADDED-HPA.xlsx" using gsub function:
fn.out = gsub(".xlsx", "_ADDED-HPA.xlsx", mytablefile); fn.out
## output the merge as a single sheet Excel file:
write.xlsx(x = df_HPA, file = fn.out)
Match example #2 (more complex because you are starting with protein names):
Fix inconsistent protein/gene names and add selected HPA columns
Suppose you are starting with a different dataframe df containing protein-protein interactions. You have columns "ProteinA" and "ProteinB" in your original dataset, but not Ensembl Gene IDs or gene names. You need to adjust for the following potential issues:
- Clarify which dataset column you are matching to the HPA dataset. Is it ProteinA or ProteinB? I recommend doing this by added a prefix or suffix to the HPA dataframe column names, e.g. "B_Gene" or "Gene_B" instead of "Gene".
- Update protein names to gene names. Most will be the same value, but not all.
## PREPARE FOR MERGE IF YOU HAVE PROTEIN NAMES ################################
## Copy hpa_subset twice and edit the column names: ------------------------
hpa_A <- hpa_subset
colnames(hpa_A) <- paste0("A_", colnames(hpa_subset))
colnames(hpa_A)
dim(hpa_A)
hpa_B <- hpa_subset
colnames(hpa_B) <- paste0("B_", colnames(hpa_subset))
colnames(hpa_B)
dim(hpa_B)
## Below I will just focus on adding matched data for ProteinB column, but you can do the same thing with ProteinA column and hpa_A dataset...
## Update your input values from ProteinB (want current gene values) -------------
## I changed gene_check to delete /// and everything afterward. This fixes an issue with
## a weird row that had Suggested.Symbol = "GPAT3 /// LPCAT1" ## check genes with R package HGNChelper
library(HGNChelper)
#gene_check <- checkGeneSymbols(df$ProteinB) ## old way (commented out)
gene_check <- checkGeneSymbols(df$ProteinB) %>%
mutate(
Suggested.Symbol = str_trim(sub("///.*", "", Suggested.Symbol)) ## error fix
)
## Check results
nrow(gene_check); nrow(df)
View(gene_check)
## Add standardized gene column for your dataset:
df <- df %>%
mutate(
ProteinB_Gene = ifelse(
gene_check$Approved, ## TRUE or FALSE values
ProteinB, ## return this value if above is TRUE
gene_check$Suggested.Symbol ## otherwise return this value
)
)
## ADD SELECTED HPA COLUMNS TO YOUR DATA #####################################df_HPA <- left_join(x = df, ## your protein-protein interaction data
y = hpa_B, ## this is hpa_subset but with column names edited
by = c("ProteinB_Gene" = "Gene"))
dim(df); dim(df_HPA); dim(hpa_B)
View(df_HPA)
## OUTPUT ####################################################################
## create output filename of your choice:
fn.out = "myoutput.xlsx"
## alternatively, create output filename from your input filename,
## but replace ".xlsx" with "_ADDED-HPA.xlsx" using gsub function:
fn.out = gsub(".xlsx", "_ADDED-HPA.xlsx", mytablefile); fn.out
## output the merge as a single sheet Excel file:
write.xlsx(x = df_HPA, file = fn.out)
Final steps and advice
Congratulations! You annotated your starting dataset with additional information from the Human Protein Atlas. This skill will save you a lot of time in future projects.
My final advice to you is this:
Always look for sources of errors. Assume the code isn't perfect.
Every time you merge data, I recommend that you "spot check" a few rows manually. Focus on missing values, any values with unusual syntax (e.g. anything that Excel might convert to a date like gene SEPT9), and values at the end. If there are any missing values, check the Human Protein Atlas website by hand and see if you can find a match without R code. Doing these manual checks is important to make sure that the code is doing what you expect.
The code in this tutorial already solved several issues that I discovered doing my own "spot checks" after merging data. For example, importing the Human Protein Atlas tab-separated value (.tsv) file in the default way like this resulted in missing rows:
hpa <- read.table("./temp/proteinatlas.tsv",
sep="\t",
header=TRUE,
stringsAsFactors=FALSE)
I opened the original proteinatlas.tsv file myself in Notepad++ to view the data as text. Notepad++ is a text viewer for Windows that works better for large files than the built-in Notepad software. Very useful for bioinformatics. I found one row with a single quotation mark that made R think that the file ended early, so my reference dataframe hpa was missing several rows. I updated my import code to solve this issue by telling R to treat quotation marks as text instead of special code characters:
hpa <- read.table("./temp/proteinatlas.tsv",
sep="\t", ## tab-separated files
quote = "\"",
header=TRUE,
stringsAsFactors=FALSE)
In another case, everything matched except a few rows, and I realized with manual checks that I had a mismatch between protein and gene symbols. Adding the extra step with R package HGNChelper helped solve that.
Always double check your work! Good luck and congratulations on the new skill.
No comments:
Post a Comment