Combined T Cell Lymphoma Cell Line Analysis of RNA-seq, PRO-seq, and ATAC-seq

.libPaths("/n/data1/cores/ntc/R-4.0.1/library/")
setwd("/n/scratch3/groups/hms/cores/ntc/Weinstock008")

library(DESeq2)

working on master data file with scripts run (bedtools closest and intersect) to link GGA annot genes to dREG filtered annot enhancers. add in counts data from PRO and ATAC across dREG span.

full <- read.table("overall.bed", sep="\t")
colnames(full) <- c("chr_gene", "start_gene", "end_gene", "name_gene", "score_gene", "strand_gene", 
     "chr_peak", "start_peak", "end_peak", "name_peak", "score_peak", "strand_peak", "biotype", "distance")

write.table(full, "TCL_Gene_to_Enhancer_master.txt", quote = FALSE, sep = "\t", row.names = FALSE, col.names = TRUE)

load DEseq and normalize with DEseq size factors based on depth Import gene body counts table generated in O2 into R

normalizing counts for all data types with DEseq size factors and melding to full table

raw_atac <- read.table("TCL_ATAC_counts.txt", header = TRUE, sep = "\t", row.names = 1)
raw_pro <- read.table("TCL_PRO_counts.txt", header = TRUE, sep = "\t", row.names = 1)
raw_rna <- read.table("TCL_RNA_counts.txt", header = TRUE, sep = "\t", row.names = 1)

ATAC first

Create metadata object to categorize samples:

metadata_atac <- data.frame(row.names = colnames(raw_atac[, 6:ncol(raw_atac)]),
                       cell_line = c(rep("TCL",23),rep("Other",1)))
                       
metadata_atac # verify the table looks ok

#confirm that row and column names match
all(rownames(metadata_atac) %in% colnames(raw_atac))
all(rownames(metadata_atac) == colnames(raw_atac)[6:ncol(raw_atac)])


## Create DESeq Data Set from the count table, and sample list:

dds_atac <- DESeqDataSetFromMatrix(countData = raw_atac[, 6:ncol(raw_atac)],      # the matrix of counts from featureCounts;
                              colData = metadata_atac,                     # the sample table created in the last step;
                              design= ~cell_line)            # this is R's notation for models

dds_atac <- estimateSizeFactors(dds_atac)
deseq_sizes_atac <- as.data.frame(sizeFactors(dds_atac))
deseq_sizes_atac
write.table(deseq_sizes_atac, "ATAC_norm_factors.txt", quote = FALSE, sep = "\t", row.names = TRUE, col.names = TRUE)

# normalization of ATAC data

dds_atac <- estimateDispersions(dds_atac)
pdf("Dispersion_ATAC_norm.pdf")
plotDispEsts(dds_atac)
dev.off()

# check results of normalization
colSums(counts(dds_atac))                # raw read counts per sample
colSums(counts(dds_atac, normalized=T))  # Total number of normalized counts per sample
data.frame(colSums(counts(dds_atac)), (colSums(counts(dds_atac, normalized=T))), row.names=colnames(counts(dds_atac)))

write.table(counts(dds_atac, normalized=TRUE), file="ATAC_norm_counts.txt", sep="\t")
counts_atac <- read.table("ATAC_norm_counts.txt", sep="\t")
counts_atac$dREG <- row.names(counts_atac)

full <- merge(full, counts_atac, by.x = "name_peak", by.y = "dREG", all.x=TRUE)

Now PRO

### Create metadata object to categorize samples:

metadata_pro <- data.frame(row.names = colnames(raw_pro[, 6:ncol(raw_pro)]),
                            cell_line = c(rep("TCL_A",12),rep("TCL_B",12)))

metadata_pro # verify the table looks ok

#confirm that row and column names match
all(rownames(metadata_pro) %in% colnames(raw_pro))
all(rownames(metadata_pro) == colnames(raw_pro)[6:ncol(raw_pro)])


## Create DESeq Data Set from the count table, and sample list:

dds_pro <- DESeqDataSetFromMatrix(countData = raw_pro[, 6:ncol(raw_pro)],      # the matrix of counts from featureCounts;
                                   colData = metadata_pro,                     # the sample table created in the last step;
                                   design= ~cell_line)            # this is R's notation for models

dds_pro <- estimateSizeFactors(dds_pro)
deseq_sizes_pro <- as.data.frame(sizeFactors(dds_pro))
deseq_sizes_pro
write.table(deseq_sizes_pro, "PRO_norm_factors.txt", quote = FALSE, sep = "\t", row.names = TRUE, col.names = TRUE)

#normalization of PRO data

dds_pro <- estimateDispersions(dds_pro)
pdf("Dispersion_PRO_norm.pdf")
plotDispEsts(dds_pro)
dev.off()

# check results of normalization
colSums(counts(dds_pro))                # raw read counts per sample
colSums(counts(dds_pro, normalized=T))  # Total number of normalized counts per sample
data.frame(colSums(counts(dds_pro)), (colSums(counts(dds_pro, normalized=T))), row.names=colnames(counts(dds_pro)))

write.table(counts(dds_pro, normalized=TRUE), file="PRO_norm_counts.txt", sep="\t")
counts_pro <- read.table("PRO_norm_counts.txt", sep="\t")
counts_pro$dREG <- row.names(counts_pro)


full <- merge(full, counts_pro, by.x = "name_peak", by.y = "dREG", all.x=TRUE)

Now RNA

### Create metadata object to categorize samples:

metadata_rna <- data.frame(row.names = colnames(raw_rna[, 6:ncol(raw_rna)]),
                            cell_line = c(rep("TCL_A",12),rep("TCL_B",12)))

metadata_rna # verify the table looks ok

#confirm that row and column names match
all(rownames(metadata_rna) %in% colnames(raw_rna))
all(rownames(metadata_rna) == colnames(raw_rna)[6:ncol(raw_rna)])


## Create DESeq Data Set from the count table, and sample list:

dds_rna <- DESeqDataSetFromMatrix(countData = raw_rna[, 6:ncol(raw_rna)],      # the matrix of counts from featureCounts;
                                   colData = metadata_rna,                     # the sample table created in the last step;
                                   design= ~cell_line)            # this is R's notation for models

dds_rna <- estimateSizeFactors(dds_rna)
deseq_sizes_rna <- as.data.frame(sizeFactors(dds_rna))
deseq_sizes_rna
write.table(deseq_sizes_rna, "RNA_norm_factors.txt", quote = FALSE, sep = "\t", row.names = TRUE, col.names = TRUE)

#normalization of RNA data

dds_rna <- estimateDispersions(dds_rna)
pdf("Dispersion_RNA_norm.pdf")
plotDispEsts(dds_rna)
dev.off()

# check results of normalization
colSums(counts(dds_rna))                # raw read counts per sample
colSums(counts(dds_rna, normalized=T))  # Total number of normalized counts per sample
data.frame(colSums(counts(dds_rna)), (colSums(counts(dds_rna, normalized=T))), row.names=colnames(counts(dds_rna)))

write.table(counts(dds_rna, normalized=TRUE), file="RNA_norm_counts.txt", sep="\t")
counts_rna <- read.table("RNA_norm_counts.txt", sep="\t")
counts_rna$dREG <- row.names(counts_rna)

full <- merge(full, counts_rna, by.x = "name_peak", by.y = "dREG", all.x=TRUE)

final table completed

write.table(full, "TCL_Final_gene_to_enhancer.txt", quote = FALSE, sep = "\t", row.names = FALSE, col.names = TRUE)

# clean-up and save out DEseq data sets

saveRDS(dds_atac, file="ATAC.rds")
saveRDS(dds_pro, file="PRO.rds")
saveRDS(dds_rna, file="RNA.rds")

######### END #########
## MKJ 24JUL24