Creative Commons License
This blog by Tommy Tang is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

My github papge

Monday, May 22, 2017

when NAs are not NAs

Using TCGAbiolinks, I downloaded RNAseq data for LUAD and LUSC

# query_rna_LUAD.hg38 <- GDCquery(project = "TCGA-LUAD", data.category = "Transcriptome Profiling",
#                   data.type = "Gene Expression Quantification", 
#                   workflow.type = "HTSeq - Counts")
# query_rna_LUSC.hg38 <- GDCquery(project = "TCGA-LUSC", data.category = "Transcriptome Profiling",
#                   data.type = "Gene Expression Quantification", 
#                   workflow.type = "HTSeq - Counts")
# GDCdownload(query_rna_LUAD.hg38, method = "client")
# GDCdownload(query_rna_LUSC.hg38, method = "client")
# LUAD_rna_data <- GDCprepare(query_rna_LUAD.hg38)
# LUSC_rna_data <- GDCprepare(query_rna_LUSC.hg38)

# I have saved both R objects into disk
# a RangedSummarizedExperiment object
## class: RangedSummarizedExperiment 
## dim: 57035 551 
## metadata(0):
## assays(1): HTSeq - Counts
## rownames(57035): ENSG00000000003 ENSG00000000005 ...
##   ENSG00000281912 ENSG00000281920
## rowData names(3): ensembl_gene_id external_gene_name
##   original_ensembl_gene_id
## colnames(551): TCGA-77-8009-01A-11R-2187-07
##   TCGA-34-5239-01A-21R-1820-07 ... TCGA-NK-A7XE-01A-12R-A405-07
##   TCGA-43-6773-11A-01R-1949-07
## colData names(69): patient barcode ...
##   subtype_Homozygous.Deletions subtype_Expression.Subtype
The problem is with the meta data
LUSC_coldata<- colData(LUSC_rna_data)
LUAD_coldata<- colData(LUAD_rna_data)

we will see the different representations of NAs

table(LUAD_coldata$subtype_Smoking.Status, useNA = "ifany")
##      Current reformed smoker for > 15 years 
##                                          78 
## Current reformed smoker for < or = 15 years 
##                                          77 
##                              Current smoker 
##                                          47 
##                         Lifelong Non-smoker 
##                                          32 
##                             [Not Available] 
##                                          11 
##                                        <NA> 
##                                         349
table(LUSC_coldata$subtype_Smoking.Status, useNA = "ifany")
##      Current reformed smoker for > 15 years 
##                                          51 
## Current reformed smoker for < or = 15 years 
##                                          87 
##                              Current smoker 
##                                          28 
##                         Lifelong Non-smoker 
##                                           7 
##                                         N/A 
##                                           6 
##                                        <NA> 
##                                         372
we see NAs are represented either as real <NA>[Not Avaiable] or N/A. The first thing to do is to tidy the metadata changing all NAs to <NA>:
I will use stringr from the tidyverse packages.
# one needs to \\ to escape [ and ]
str_replace(LUAD_coldata$subtype_Smoking.Status, "\\[Not Available\\]", NA)
## Error: `replacement` must be a character vector
This does not quite work. stringr enforces that the replacement is the same type as a character. instead:
str_replace(LUAD_coldata$subtype_Smoking.Status, "\\[Not Available\\]", NA_character_) %>% table(useNA = "ifany")
## .
## Current reformed smoker for < or = 15 years 
##                                          77 
##      Current reformed smoker for > 15 years 
##                                          78 
##                              Current smoker 
##                                          47 
##                         Lifelong Non-smoker 
##                                          32 
##                                        <NA> 
##                                         360
str_replace(LUSC_coldata$subtype_Smoking.Status, "N/A", NA_character_) %>% table(useNA = "ifany")
## .
## Current reformed smoker for < or = 15 years 
##                                          87 
##      Current reformed smoker for > 15 years 
##                                          51 
##                              Current smoker 
##                                          28 
##                         Lifelong Non-smoker 
##                                           7 
##                                        <NA> 
##                                         378
to test if a vector contains any NAs:
# this is invalid
# LUAD_coldata$subtype_Smoking.Status == NA, instead, use$subtype_Smoking.Status) %>% table()
## .
##   245   349
The traditional way to change a string to NA is:
myvector<- c("NA", "NA", "a", "b", "c", NA)
## [1] "NA" "NA" "a"  "b"  "c"  NA
myvector[myvector =="NA"]<- NA
## [1] NA  NA  "a" "b" "c" NA


  1. Data are messy, even if the data are packaged from a big project such as TCGA.
  2. Different representation of NA could be a devil in your data, tidy it to R’s representation is needed.
  3. stringr is a great tool in your belt for tidying data.

Thursday, April 20, 2017

SHERLOCK to detect early circulating tumor DNA

Early detection of circulating tumor DNA is quite challenging. Read this comment in Cell:

And here comes SHERLOCK(Specific High Sensitivity Enzymatic Reporter UnLOCKing)!
from Nucleic acid detection with CRISPR-Cas13a/C2c2,  a study leaded by Feng Zhang in MIT.

Rapid, inexpensive, and sensitive nucleic acid detection may aid point-of-care pathogen detection, genotyping, and disease monitoring. The RNA-guided, RNA-targeting CRISPR effector Cas13a (previously known as C2c2) exhibits a “collateral effect” of promiscuous RNAse activity upon target recognition. We combine the collateral effect of Cas13a with isothermal amplification to establish a CRISPR-based diagnostic (CRISPR-Dx), providing rapid DNA or RNA detection with attomolar sensitivity and single-base mismatch specificity. We use this Cas13a-based molecular detection platform, termed SHERLOCK (Specific High Sensitivity Enzymatic Reporter UnLOCKing), to detect specific strains of Zika and Dengue virus, distinguish pathogenic bacteria, genotype human DNA, and identify cell-free tumor DNA mutations. Furthermore, SHERLOCK reaction reagents can be lyophilized for cold-chain independence and long-term storage, and readily reconstituted on paper for field applications.

This may be a game-changer in the field.

several very interesting histone modification (epigenetics) papers

If you have read my previous post Misconcept of epigenetics , you will know that: histone modification is not epigenetics!

However, recently, four papers in Science actually showed that histone modifications actually have a component of inheritance:

Causal role for inheritance of H3K27me3 in maintaining the OFF state of a Drosophila HOX gene

DNA sequence-dependent epigenetic inheritance of gene silencing and histone H3K9 methylation

Propagation of Polycomb-repressed chromatin requires sequence-specific recruitment to DNA

Transgenerational transmission of environmental information in C. elegans

and one paper in Molecular Cell:

SNF2 Family Protein Fft3 Suppresses Nucleosome Turnover to Promote Epigenetic Inheritance and Proper Replication

one more on Nature Genetics:

Interestingly, all the studies were carried out on model organisms: Drosophila and C.elegans!

Read this blog post if you are into this question: what is epigenetics ?

Wednesday, April 19, 2017

sshfs on ubuntu and ssh key

Two things I want to keep a note here:

First, if you ever have set up a shh key for connecting to remote server, you need to be aware that password-less shh key only works when your home directory on the server is not 777 (writable by others).

see this stackexchange post .

Second, I was following to set up sshfs on my ubuntu machine. I put down a gist below.

Friday, February 24, 2017

use bash associate array to generate commands

The problem

I am running pyclone recently to determine the clonality of our collected tumor samples. It needs tumor purity (estimated from sequenza) for input. I have 24 samples, and I want to generate pyclone commands for all of them.

The solution

I usually generate command by bash and then use another bash wrapper to generate pbs files on HPC ((Thanks to @SBAmin). now for each patient, I have two samples. How should I generate the commands? I am still better in R and Unix than python, so I used the associated array in bash.
First, generate a file containing the tumor purity for each tumor:
head -6 all_tumor_purity_no_header.txt
0.69    Pa25T1
0.26    Pa25T2
0.49    Pa26T1
0.37    Pa26T2
0.9 Pa27T1
0.92    Pa27T2
This bash script uses associate array to contain tumor purity. read more at
#! /bin/bash
set -euo pipefail

## build the array to contain the tumor purity, like a python dict
## have to declare by -A
declare -A cols

while read purity sample
done < all_purity_no_header.txt 

echo ${cols[@]}

## generate commands
for i in Pa{25..37}
   echo PyClone run_analysis_pipeline --in_file ${i}T1_pyclone.tsv ${i}T2_pyclone.tsv --tumour_contents ${cols[${i}T1]} ${cols[${i}T2]} --samples ${i}T1 ${i}T2 --density pyclone_binomial --working_dir ${i}T_pyclone_analysis --min_cluster_size 2 --seed 123 --num_iters 50000 > ${i}_pyclone_commands.txt
chmod u+x
what you get: cat *commands.txt
PyClone run_analysis_pipeline --in_file Pa25T1_pyclone.tsv Pa25T2_pyclone.tsv --tumour_contents 0.69 0.26 --samples Pa25T1 Pa25T2 --density pyclone_binomial --working_dir Pa25T_pyclone_analysis --min_cluster_size 2 --seed 123 --num_iters 50000
PyClone run_analysis_pipeline --in_file Pa26T1_pyclone.tsv Pa26T2_pyclone.tsv --tumour_contents 0.49 0.37 --samples Pa26T1 Pa26T2 --density pyclone_binomial --working_dir Pa26T_pyclone_analysis --min_cluster_size 2 --seed 123 --num_iters 50000
PyClone run_analysis_pipeline --in_file Pa27T1_pyclone.tsv Pa27T2_pyclone.tsv --tumour_contents 0.9 0.92 --samples Pa27T1 Pa27T2 --density pyclone_binomial --working_dir Pa27T_pyclone_analysis --min_cluster_size 2 --seed 123 --num_iters 50000
use the makemsub wrapper:

find *commands.txt | parallel 'makemsub -a {} -j {/.} -o a -c 4 -t "48:00:00" -m 16g > {/.}.pbs'

for pbs in *pbs
  msub $pbs
  sleep 1
I usually do it for simple tasks. e.g. only one or two commands are evoked. For more complex workflows, a workflow tool such as snakemake is better.

Wednesday, February 1, 2017

How to enable the scroll mode for tmux

tmux config file

it changed the key binding from control + b to control + a if you are familiar with the screen shortcuts.
control + a + c will create a new window. control + a + Space will move to previous window. control + a + n will move to next window.
control + a + ?
will show you all the shortcuts.

scroll mode

One problem with screen or tmux is that you have to press control + a + [ to enter the copy mode, and and control + a + ] to paste it. I want to just use the mouse to scroll up and down and copy/paste.
read this long thread github issue: The solution that worked for me:
git clone ~/.tmux/plugins/tpm
Put this at the bottom of .tmux.conf:
# List of plugins
set -g @plugin 'tmux-plugins/tpm'
set -g @plugin 'tmux-plugins/tmux-sensible'

# Other examples:
# set -g @plugin 'github_username/plugin_name'
# set -g @plugin ''
# set -g @plugin ''

# Initialize TMUX plugin manager (keep this line at the very bottom of tmux.conf)
run '~/.tmux/plugins/tpm/tpm'

open your .tmux.conf.
To enable mouse-mode in tmux 2.1+, put the following line in your ~/.tmux.conf:
set-option -g mouse on

then add the following line to your .tmux.conf file:
set -g @plugin 'nhdaly/tmux-better-mouse-mode'
  • install it
# start a new session

# install plugin
`control + a + I (captial)` to install all the plugins.

Now if you scroll up with your mouse, you will enter into copy mode automatically, and when you scroll down to the end of the current screen, you will exit the copy mode automatically.
  • copy and paste
If you scroll up and select the text you want to copy by left-click and drag, you will exit the copy mode instantly, and the content you selected will be copied in the buffer. You just need to control + a + ] to paste it. very cool!

Monday, January 30, 2017

How to make a heatmap using character matrix

People are willing to read colors. I will show you how to repesent characters in heatmap. ComplexHeatmap can directly translate a character matrix into heatmap. It is good to learn basics of grid package which ComplexHeatmap and ggplot2 build on.
grid.rect is the function used to draw rectangles.
colors<- structure(circlize::rand_color(4), names = c("a", "b", "c", "d"))
mat_letters<- matrix(sample(letters[1:4], 100, replace = TRUE), 10, 10)
## default
Heatmap(mat_letters, col = colors)
I want to have some gap between the cells:
cell_fun<-  function(j, i, x, y, width, height, fill) {
    grid.rect(x = x *0.6, y = y, width = width * 0.6, height = height, 
        gp = gpar(col = "grey", fill = fill))
Heatmap(mat_letters, rect_gp = gpar(type = "none"), cell_fun = cell_fun, col = colors)
It may look different in an interactive graphical device. save the figure in a pdf and you will see the effect.
## always square and have gaps between cells
cell_fun2<- function(j, i, x, y, width, height, fill) {
    s = min(unit.c(convertWidth(width, "cm"), convertHeight(height, "cm")))
    grid.rect(x = x*0.8, y = y*0.8, width = s * 0.8, height = s*0.8, 
        gp = gpar(col = "grey", fill = fill))
Heatmap(mat_letters, rect_gp = gpar(type = "none"), cell_fun = cell_fun2, col = colors)
cell_fun3<-  function(j, i, x, y, width, height, fill) {
    grid.rect(x = x*0.8, y = y, width = width * 0.8, height = height, 
        gp = gpar(col = "grey", fill = fill))
Heatmap(mat_letters, rect_gp = gpar(type = "none"), cell_fun = cell_fun3, col = colors)

make cells thiner

cell_fun4<- function(j, i, x, y, width, height, fill) {
        s = min(unit.c(convertWidth(width, "cm"), convertHeight(height, "cm")))
    grid.rect(x = x*0.6 , y = y, width = s* 0.6, height = height, 
        gp = gpar(col = "grey", fill = fill))
Heatmap(mat_letters, rect_gp = gpar(type = "none"), cell_fun = cell_fun4, col = colors, show_heatmap_legend = FALSE)