-
Notifications
You must be signed in to change notification settings - Fork 1
/
08_tximport_aggregate_transcript_counts.R
100 lines (65 loc) · 2.9 KB
/
08_tximport_aggregate_transcript_counts.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
#####################################################################################################################################################################
#
# Script used in:
#
# "Gene expression profiling across ontogenetic stages in wood white (Leptidea sinapis) reveals pathways linked to butterfly diapause regulation."
# Luis Leal, Venkat Talla, Thomas Källman, Magne Friberg, Christer Wiklund, Vlad Dincă, Roger Vila, Niclas Backström
# Mol. Eco. (2018)
#
#
#####################################################################################################################################################################
#
# Luis Leal
# Uppsala University, Uppsala, Sweden, 2017
#
#####################################################################################################################################################################
###################################################################################################
#
# R/3.4.0 script used to convert Kallisto output files from count-per-transcript to count-per-gene
#
# Dependencies: abundance.tsv >> kallisto output file
# samples_ID_tximport.txt >> basic identification of each sample (experimental design)
# gene_names.csv >> associates each transcript to its parental gene
#
###################################################################################################
###Install tximport
#source("https://bioconductor.org/biocLite.R")
#biocLite("tximport")
### Install HDF5 interface to R package
#biocLite("rhdf5")
#### Install readr package
#install.packages("readr")
###load tximport and rhdf5 libraries
library(tximport)
library(rhdf5)
library(readr)
###Clear all states
rm(list=ls(all=TRUE))
##### PATHS
### locate the directory containing the kallisto files; list files in directory
base_dir <- "/kallisto_OUT"
list.files(base_dir)
### Specify where auxiliary table describing the experimental design is stored
table_dir <- "/tximport"
### Specify name of auxiliary table describing the experimental design
auxTable_filename <- "samples_ID_tximport.txt"
### Specify output folder
out_dir <- "/tximportOUT"
##### Import auxiliary files
### Read auxiliary table; show table
samples <- read.table(file.path(table_dir, auxTable_filename), header = TRUE)
samples
### Read file that associates transcripts with gene IDs (necessary for gene-level summarization)
tx2gene <- read.csv(file.path(table_dir, "gene_names.csv"), header = TRUE)
head(tx2gene)
##### Import Kallisto file
files <- file.path(base_dir, "results", samples$experiment, "abundance.tsv")
names(files) <-t(samples[1])
files
##### Run tximport
### Note: do NOT normalize by library size (DESeq2 will perform this step)
txi <- tximport(files, type = "kallisto", tx2gene = tx2gene)
names(txi)
head(txi$counts)
# Write data to file in CSV format
write.table(txi$counts, file.path(out_dir, "tximport_OUTPUT.csv"), sep="\t")