EM_test/RSEM
parallel -j 2 rsem-calculate-expression --bowtie2 -p 4 {} ../Reference_Genome/hg19 {/.} ::: ../FASTQ/*.fastq
```

We can now put all the RSEM counts/TPM results in one matrix, if it hasn't been done yet.


```r
rsem_files <- list.files("RSEM_test/RSEM/", pattern="genes.results", full.names = TRUE)

# Read all files and create a list of data.tables
rsem_list <- lapply(rsem_files, function(x, ...){
    y<-fread(x, drop=c("length", "effective_length", "FPKM")); 
    y$sample_name <- gsub(".genes.results", "", basename(x));
    return(y);})

# Bind all the data.tables to create a long data.table
rsem_data_long <- rbindlist(rsem_list)

# Rename a column, so that it's a bit R friendly
setnames(rsem_data_long, "transcript_id(s)", "transcript_ids")

# Reshape the long data.table to create a matrix
# Assume that missing entries would have a TPM value of 0
rsem_tpm_matrix <- dcast.data.table(rsem_data_long, gene_id+transcript_ids~sample_name, value.var = "TPM", fill=0)
rsem_count_matrix <- dcast.data.table(rsem_data_long, gene_id+transcript_ids~sample_name, value.var = "expected_count", fill=0)

# Keep track of the gene_id - transcript mapping 
rsem_txs_table <- rsem_tpm_matrix[, c("gene_id","transcript_ids"), with=FALSE]
#save(rsem_txs_table, file="RSEM_test/RSEM/rsem_txs_table.Rda")

# Remove the transcript column (we create an annotation table in the next chunk)
rsem_tpm_matrix <- rsem_tpm_matrix[, transcript_ids:=NULL][order(gene_id)]
rsem_count_matrix <- rsem_count_matrix[, transcript_ids:=NULL][order(gene_id)]

# Write to disc
write.csv(rsem_tpm_matrix, file="RSEM_test/RSEM/rsem_tpm_matrix.csv", row.names=FALSE)
write.csv(rsem_count_matrix, file="RSEM_test/RSEM/rsem_count_matrix.csv", row.names=FALSE)
```

### Get sample metadata


```r
# Write the pdata data
write.csv(pd_small, file="RSEM_test/RSEM/rsem_pdata.csv", row.names=FALSE)
```


## Get transcript annotation

Get feature annotations from BioC and create a feature annotation set file.


```r
# Install TxDb packages if you don't already have it.
t