setwd('/Volumes/BB_Backup_3//BD_aging_project/RNAseq/All_tissues_analysis/TE_analysis/TEtranscript//DE_analysis//TEST_all_genes_TE_results')
source('TE_RNAseq_analysis_functions.R')

# 2017-01-04:  RNAseq analysis for TE
#use count tables generated by TE_transcript

# 2017-03-15: 
# rerun using the new count tables generated from STAR with 150 multimapping alignments
my.TE.names <- read.table('TE_names_list.txt',header=F,sep="\t",stringsAsFactors=F)

# 2017-09-12
# new displays for paper

####################################    Liver    #################################### 
# read in TEtranscript count matrices (only 2 conditions can be measured at one time)
my.liver1 <- read.table('TE-transcript_output/Liver_29vs3m_TEtranscripts.cntTable',header=T,sep="\t",stringsAsFactors=F)
my.liver2 <- read.table('TE-transcript_output/Liver_12vs3m_TEtranscripts.cntTable',header=T,sep="\t",stringsAsFactors=F)

# verify the counts and gene names are consistent
####sum(my.liver1$gene.TE %in% my.liver2$gene.TE) # 23514
####sum(my.liver1$..NEW_STAR_map.3m1_Liver_ATCACG_RAligned.out.bam.C  %in% my.liver2$..NEW_STAR_map.3m1_Liver_ATCACG_RAligned.out.bam.C) # 23514

my.liver <- cbind(my.liver2[,c(1,5:7,2:4)],my.liver1[,2:4])
rownames(my.liver) <- my.liver[,1]

# process RNAseq data and save RData object
my.liver.RNAseq.process <- process_TE_rnaseq("Liver", my.liver,my.TE.names = my.TE.names)
save(my.liver.RNAseq.process, file=paste(Sys.Date(),"RNA_seq_result_Liver.RData", sep="_") )
##################################################################################### 

####################################   Heart   #################################### 
# read in TEtranscript count matrices
my.heart1 <- read.table('TE-transcript_output/Heart_29vs3m_TEtranscripts.cntTable',header=T,sep="\t",stringsAsFactors=F)
my.heart2 <- read.table('TE-transcript_output/Heart_12vs3m_TEtranscripts.cntTable',header=T,sep="\t",stringsAsFactors=F)
my.heart <- cbind(my.heart2[,c(1,5:7,2:4)],my.heart1[,2:4])
rownames(my.heart) <- my.heart[,1]

# process RNAseq data and save RData object
my.heart.RNAseq.process <- process_TE_rnaseq("Heart", my.heart,my.TE.names = my.TE.names)
save(my.heart.RNAseq.process, file=paste(Sys.Date(),"RNA_seq_result_Heart.RData", sep="_") )
###################################################################################

#################################### Cerebellum #################################### 
# read in TEtranscript count matrices
my.cereb1 <- read.table('TE-transcript_output/Cerebellum_29vs3m_TEtranscripts.cntTable',header=T,sep="\t",stringsAsFactors=F)
my.cereb2 <- read.table('TE-transcript_output/Cerebellum_12vs3m_TEtranscripts.cntTable',header=T,sep="\t",stringsAsFactors=F)

# collapse the 2 runs
my.cereb <- data.frame(gene.TE = my.cereb1$gene.TE,
                       '3m1'= (my.cereb1$..NEW_STAR_map.Cereb_3m1_1st_run_ATCACG_RAligned.out.bam.C + my.cereb1$..NEW_STAR_map.Cereb_3m1_2nd_run_ATCACG_RAligned.out.bam.C),
                       '3m2' = (my.cereb1$..NEW_STAR_map.Cereb_3m2_1st_run_CGATGT_RAligned.out.bam.C + my.cereb1$..NEW_STAR_map.Cereb_3m2_2nd_run_CGATGT_RAligned.out.bam.C),
                       '3m3' =  (my.cereb1$..NEW_STAR_map.Cereb_3m3_1st_run_TTAGGC_RAligned.out.bam.C + my.cereb1$..NEW_STAR_map.Cereb_3m3_2nd_run_TTAGGC_RAligned.out.bam.C),
                       '12m1' =  (my.cereb2$..NEW_STAR_map.Cereb_12m1_1st_run_TGACCA_RAligned.out.bam.T + my.cereb2$..NEW_STAR_map.Cereb_12m1_2nd_run_TGACCA_RAligned.out.bam.T),
                       '12m2' =  (my.cereb2$..NEW_STAR_map.Cereb_12m2_1st_run_ACAGTG_RAligned.out.bam.T + my.cereb2$..NEW_STAR_map.Cereb_12m2_2nd_run_ACAGTG_RAligned.out.bam.T),
                       '12m3' =  (my.cereb2$..NEW_STAR_map.Cereb_12m3_1st_run_GCCAAT_RAligned.out.bam.T + my.cereb2$..NEW_STAR_map.Cereb_12m3_2nd_run_GCCAAT_RAligned.out.bam.T),
                       '29m2' =  (my.cereb1$..NEW_STAR_map.Cereb_29m2_1st_run_CAGATC_RAligned.out.bam.T + my.cereb1$..NEW_STAR_map.Cereb_29m2_2nd_run_CAGATC_RAligned.out.bam.T),
                       '29m3' =   (my.cereb1$..NEW_STAR_map.Cereb_29m3_1st_run_ACTTGA_RAligned.out.bam.T + my.cereb1$..NEW_STAR_map.Cereb_29m3_2nd_run_ACTTGA_RAligned.out.bam.T),
                       '29m4' =   (my.cereb1$..NEW_STAR_map.Cereb_29m4_1st_run_GATCAG_RAligned.out.bam.T + my.cereb1$..NEW_STAR_map.Cereb_29m4_2nd_run_GATCAG_RAligned.out.bam.T)
)
rownames(my.cereb) <- my.cereb[,1]


# process RNAseq data and save RData object
my.cereb.RNAseq.process <- process_TE_rnaseq("Cerebellum", my.cereb,my.TE.names = my.TE.names)
save(my.cereb.RNAseq.process, file=paste(Sys.Date(),"RNA_seq_result_Cerebellum.RData", sep="_") )
####################################################################################

#################################### Olfactory Bulb #################################
# one of the 12mths samples was not analyzed
# read in TEtranscript count matrices
my.ob1 <- read.table('TE-transcript_output/OB_29vs3m_TEtranscripts.cntTable',header=T,sep="\t",stringsAsFactors=F)
my.ob2 <- read.table('TE-transcript_output/OB_12vs3m_TEtranscripts.cntTable',header=T,sep="\t",stringsAsFactors=F)
my.ob <- cbind(my.ob2[,c(1,4:6,2:3)],my.ob1[,2:4])
rownames(my.ob) <- my.ob[,1]

# process RNAseq data and save RData object
my.ob.RNAseq.process <- process_TE_rnaseq("OlfactoryBulb", my.ob, reps.3=3, reps.12=2, reps.29=3,my.TE.names = my.TE.names)
save(my.ob.RNAseq.process, file=paste(Sys.Date(),"RNA_seq_result_OB.RData", sep="_") )
##################################################################################### 

####################################  NPCs pools  ###################################
# read in TEtranscript count matrices
my.npc1 <- read.table('TE-transcript_output/NPCs_29vs3m_TEtranscripts.cntTable',header=T,sep="\t",stringsAsFactors=F)
my.npc2 <- read.table('TE-transcript_output/NPCs_12vs3m_TEtranscripts.cntTable',header=T,sep="\t",stringsAsFactors=F)
my.npc <- cbind(my.npc2[,c(1,4:5,2:3)],my.npc1[,2:3])
rownames(my.npc) <- my.npc[,1]

# process RNAseq data and save RData object
my.npc.RNAseq.process <- process_TE_rnaseq("NPCs", my.npc, reps.3=2, reps.12=2, reps.29=2,my.TE.names = my.TE.names)
save(my.npc.RNAseq.process, file=paste(Sys.Date(),"RNA_seq_result_NPCs.RData", sep="_") )
##################################################################################### 

