pre_process

0.1 Parameters

suffix = ""
data_to_read = ""

0.2 functions

0.3 Data

1 Pre-process

#save before omitting
genes_names = counts$gene_name
genes_lengths = counts[,5]

#omit non relevant cols
counts=counts[,7:ncol(counts)]


#create genes names
genes_names=make.unique(genes_names) %>% replace_na('NA')
rownames(counts) = genes_names

#omit non relevant genes
omitgenes= startsWith(rownames(counts),"NA")
counts=counts[!omitgenes,]
genes_lengths = genes_lengths[!omitgenes] #update genes_lengths

#calcualte gene length and MT genes
mt_genes = startsWith(rownames(counts),"MT-")| startsWith(rownames(counts),"ERCC-")

#get colnames
cell.labels <- gsub(pattern = ".bam",replacement = "",colnames(counts))

#change colnames
colnames(counts) <- cell.labels


counts_seurat <- CreateSeuratObject(counts = counts, project = "egfr_counts", min.cells = 0, min.features = 1000)

cell_type = str_extract(cell.labels, "^[A-Z]{1,3}[0-9]{3}")
condition = str_extract(cell.labels, "cp|osi|ctrl|roxa|op")
replicate = str_extract(cell.labels, ".$")
metadata = data.frame(cell_type = cell_type, condition = condition, replicate = replicate, row.names = colnames(counts))

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = round(counts),
                              colData = metadata,
                              design = ~condition)

nrow(dds)

[1] 60445

dds1 <- dds[ rowSums(counts(dds)) >= 3, ]
nrow(dds1)

[1] 36882

vst = vst(dds1, blind=FALSE)

library("ggfortify")
PCAdata <- prcomp(t(assay(vst)))
autoplot(PCAdata, data = metadata,colour = "condition",label = FALSE, main="PCA") # Show dots

Warning in grSoftVersion() :
  unable to load shared object '/usr/local/lib/R/modules//R_X11.so':
  libXt.so.6: cannot open shared object file: No such file or directory

VlnPlot(counts_seurat, features = c("nCount_RNA"),group.by = "orig.ident",pt.size = 3)+theme(text = element_text(size=10),
        axis.text.x = element_text(size = 8)) 

Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Computation failed in `stat_ydensity()`:
replacement has 1 row, data has 0

counts_filtered = counts[, colnames(counts) != "HCC827cp3"]
metadata_filtered = metadata[ rownames(metadata) != "HCC827cp3",]
dds <- DESeqDataSetFromMatrix(countData = round(counts_filtered),
                              colData = metadata_filtered,
                              design = ~condition)

converting counts to integer mode
Warning in DESeqDataSet(se, design = design, ignoreRank) :
  some variables in design formula are characters, converting to factors

vst = vst(dds1, blind=FALSE)

PCAdata <- prcomp(t(assay(vst)))
autoplot(PCAdata, data = metadata,colour = "condition",label = FALSE, main="PCA")

autoplot(PCAdata, data = metadata,colour = "cell_type",label = FALSE, main="PCA")

pheatmap(sampleDistMatrix,
         clustering_distance_rows=sampleDists,
         clustering_distance_cols=sampleDists,
         col=colors, main = "Sample Distance Matrix ",show_colnames = T)+
    theme(axis.text.x = element_text(angle=30, hjust=1)) 

NULL