extract all genes by giving GO-terms

# # Using Bioconductor in R
# library(org.Mm.eg.db)
# genes <- select(org.Mm.eg.db, keys="GO:0048813", columns=c("SYMBOL", "ENSEMBL"), keytype="GO")
# > genes

extract annotation.gtf for mouse.

• (Optional) UCSC Genome Browser: Navigate to the UCSC Table Browser. Choose the following settings: group: “Genes and Gene Predictions” track: “RefSeq Genes” (or another gene prediction track if you have a preference) output format: “GTF – gene transfer format” Click “get output”, and you can then save the result as a .gtf file.

• (Optional) Ensembl: Go to the Ensembl FTP site. Navigate through: release-XYZ (choose the latest release) > gtf > mus_musculus (choose the appropriate assembly if there are multiple). Download the GTF file, which might be named something like Mus_musculus.GRCm38.XYZ.gtf.gz. Make sure to unzip the file after downloading. GENCODE:

• (Used) Visit the GENCODE Mouse website: https://www.gencodegenes.org/mouse/release_M25.html Download the comprehensive gene annotation GTF file for the mm10 assembly.

• Convert the GTF file to a BED file with only the TSS positions.

  awk ‘OFS=”\t” {if ($3 == “gene”) print $1, $4, $4+1, $10, $7, $12, $14}’ gencode.vM25.annotation.gtf > tss_.bed #55401 (Note not under kate)

  grep “proteincoding” tss.bed > tss.bed #21859 records

In many genomics studies, the promoter region of a gene is typically considered to be the sequence immediately upstream of the transcription start site (TSS). However, the exact definition of the “promoter region” can vary between studies and depends on the specific context. For the mouse genome (mm10), a commonly used range for the promoter region is:

• (Used) -2000 to +500 bp relative to the TSS: This defines a 2.5 kb region centered around the TSS. The 2 kb upstream often captures many of the cis-regulatory elements, while the 500 bp downstream includes the TSS and possibly the beginning of the gene. Some studies might define a narrower region, such as:

• (Optional) -1000 to +200 bp relative to the TSS Or even a broader one:

• (Optional) -5000 to +1000 bp relative to the TSS The best definition often depends on the specific research question. If the goal is to capture as many potential regulatory elements as possible, a broader range might be chosen. If the goal is to focus specifically on elements directly at the TSS, a narrower range might be more appropriate.

Convert TSS positions to promoter regions (Note that $7 is SYMBOL, $4 ENSEMBL-ID).

# For genes on the '+' strand:
awk 'OFS="\t" {start=$2-2000>0?$2-2000:0; end=$2+500; if ($5 == "+") print $1, start, end, $4}' tss.bed > promoters_plus.bed
# For genes on the '-' strand:
awk 'OFS="\t" {start=$3-500>0?$3-500:0; end=$3+2000; if ($5 == "-") print $1, start, end, $4}' tss.bed > promoters_minus.bed
# Combine the two promoter files.
cat promoters_plus.bed promoters_minus.bed > promoters.bed

# #(Optional) Extract sequences for these promoter regions.
# bedtools getfasta -fi mm10_genome.fasta -bed promoters.bed -fo promoters.fasta

Quantify Methylation Signal in Promoters: With the promoter regions defined, use your bam files to calculate coverage or read depth over these regions. The depth can be interpreted as a signal of methylation. You can use tools such as bedtools coverage to get the depth of coverage over these promoter regions.

#Control females:
bedtools coverage -a promoters.bed -b ../alg/1_0_1_sorted.bam > promoters_coverage_1.0.1.txt
bedtools coverage -a promoters.bed -b ../alg/1_0_2_sorted.bam > promoters_coverage_1.0.2.txt
bedtools coverage -a promoters.bed -b ../alg/1_0_3_sorted.bam > promoters_coverage_1.0.3.txt
#Control males:
bedtools coverage -a promoters.bed -b ../alg/1_1_2b_sorted.bam > promoters_coverage_1.1.2b.txt
bedtools coverage -a promoters.bed -b ../alg/1_1_3_sorted.bam > promoters_coverage_1.1.3.txt
bedtools coverage -a promoters.bed -b ../alg/1_1_4_sorted.bam > promoters_coverage_1.1.4.txt
#Stress females:
bedtools coverage -a promoters.bed -b ../alg/2_0_1b_sorted.bam > promoters_coverage_2.0.1b.txt
bedtools coverage -a promoters.bed -b ../alg/2_0_2_sorted.bam > promoters_coverage_2.0.2.txt
bedtools coverage -a promoters.bed -b ../alg/2_0_3_sorted.bam > promoters_coverage_2.0.3.txt
#Stress males:
bedtools coverage -a promoters.bed -b ../alg/2_1_1_sorted.bam > promoters_coverage_2.1.1.txt
bedtools coverage -a promoters.bed -b ../alg/2_1_2_sorted.bam > promoters_coverage_2.1.2.txt
bedtools coverage -a promoters.bed -b ../alg/2_1_3_sorted.bam > promoters_coverage_2.1.3.txt

Parse and consolidate coverage files:

# delete the 'chr positions '
for file in ./promoters_coverage_*.txt; do
    awk '{$1=$2=$3=""; print substr($0, 4)}' $file > tmp.txt && mv tmp.txt $file
done
for file in ./promoters_coverage_*.txt; do
    sed 's/^ *//' $file > tmp.txt && mv tmp.txt $file
done
# text if they promoters_coverage_*.txt contain exactly the same set of genes (and in the same order),
cut -d " " -f1 promoters_coverage_1.0.1.txt > reference_genes.txt
for file in ./promoters_coverage_*.txt; do
    if ! cmp -s <(cut -d " " -f1 $file) reference_genes.txt; then
        echo "$file does not match reference!"
    fi
done
#merge the counts
cut -d " " -f1 promoters_coverage_1.0.1.txt > read_counts_table.txt
for file in ./promoters_coverage_*.txt; do
    paste read_counts_table.txt <(cut -d " " -f2 $file) > tmp && mv tmp read_counts_table.txt
done
#add the sample names
for file in ./promoters_coverage_*.txt; do
    basename $file | sed 's/promoters_coverage_//;s/.txt//'
done > sample_names.txt
echo -n "Gene " > header.txt
paste -s -d ' ' sample_names.txt >> header.txt
{
    cat header.txt
    cat read_counts_table.txt
} > read_counts_table_with_header.txt

#MANUALLY replace "\n" in the first line, DELETE " and "; 
#DELETE the version of EnsembleID with the following command: 21860-1 records remaining.
cp read_counts_table_with_header.txt temp
awk 'BEGIN {FS=OFS="\t"} NR==1 {print; next} {gsub(/\.+[0-9]+$/, "", $1); print}' temp > read_counts_table_with_header.txt

library(DESeq2)
setwd("/home/jhuang/DATA/Data_Emilia_MeDIP/analysis_2023_3")
# Load the data
count_data <- read.table("read_counts_table_with_header.txt", header = TRUE, sep = "\t", quote = "", row.names = 1)
sampleInfo <- read.csv("sampleInfo.txt", header = TRUE, stringsAsFactors = FALSE)
row.names(sampleInfo) <- sampleInfo$sampleID
dds <- DESeqDataSetFromMatrix(
  countData = count_data,
  colData = sampleInfo,
  design = ~ condition
)

# Plot PCA
dds <- DESeq(dds)
rld <- rlogTransformation(dds)
png("pca.png", 1200, 800)
plotPCA(rld, intgroup=c("condition"))
dev.off()

# Differential analysis
dds$condition <- relevel(dds$condition, "control_female")
dds = DESeq(dds, betaPrior=FALSE)
resultsNames(dds)
clist <- c("stress_female_vs_control_female")

dds$condition <- relevel(dds$condition, "control_male")
dds = DESeq(dds, betaPrior=FALSE)
resultsNames(dds)
clist <- c("stress_male_vs_control_male")

library(biomaRt)
listEnsembl()
listMarts()
listDatasets(ensembl)
ensembl <- useEnsembl(biomart = "ensembl", dataset = "mmusculus_gene_ensembl", version="100")
datasets <- listDatasets(ensembl)
listEnsemblArchives()
attributes = listAttributes(ensembl)
attributes[1:25,]
for (i in clist) {
  #i<-clist[1]
  contrast = paste("condition", i, sep="_")
  res = results(dds, name=contrast)
  res <- res[!is.na(res$log2FoldChange),]
  gene_ids_no_version <- gsub("\\.\\d+$", "", rownames(res))
  geness <- getBM(attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype', 'entrezgene_id', 'chromosome_name', 'start_position', 'end_position', 'strand', 'description'),
      filters = 'ensembl_gene_id',
      values = gene_ids_no_version, 
      mart = ensembl)
  geness_uniq <- distinct(geness, ensembl_gene_id, .keep_all= TRUE)
  #merge by column by common colunmn name, in the case "GENEID"
  res$ENSEMBL = rownames(res)
  identical(rownames(res), geness_uniq$ensembl_gene_id)
  res_df <- as.data.frame(res)
  geness_res <- merge(geness_uniq, res_df, by.x="ensembl_gene_id", by.y="ENSEMBL")
  dim(geness_res)
  rownames(geness_res) <- geness_res$ensembl_gene_id
  geness_res$ensembl_gene_id <- NULL
  write.csv(as.data.frame(geness_res[order(geness_res$pvalue),]), file = paste(i, "all.txt", sep="-"))
  up <- subset(geness_res, padj<=0.05 & log2FoldChange>=2)
  down <- subset(geness_res, padj<=0.05 & log2FoldChange<=-2)
  write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up.txt", sep="-"))
  write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down.txt", sep="-"))
}
#~/Tools/csv2xls-0.4/csv_to_xls.py stress_female_vs_control_female-all.txt stress_female_vs_control_female-up.txt stress_female_vs_control_female-down.txt -d',' -o stress_female_vs_control_female.xls
#~/Tools/csv2xls-0.4/csv_to_xls.py stress_male_vs_control_male-all.txt stress_male_vs_control_male-up.txt stress_male_vs_control_male-down.txt -d',' -o stress_male_vs_control_male.xls

# Clustering the genes and draw heatmap for DEGs (limiting the genes with ENSEMBL by giving GO-term or give all genes via the file 'ids_')
#install.packages("gplots")
library("gplots")

#(Optional) cat *.id | sort -u > ids  #add Gene_Id in the first line, delete the ""
#(Used) cut -d$'\t' -f1 read_counts_table_with_header.txt > ids  #Rename the first line to Gene_Id

GOI <- read.csv("ids")$Gene_Id  #21859 genes
RNASeq.NoCellLine <- assay(rld)
# Defining the custom order
#column_order <- c(
#  "uninfected.2h_r1", "uninfected.2h_r2", "uninfected.3h_r3", 
#  "uninfected.3d_r1", "uninfected.3d_r2", 
#  "1457-M10.2h_r1", "1457-M10.2h_r2", "1457-M10.2h_r3","1457.2h_r1", "1457.2h_r2", "1457.2h_r3",  
#  "1457-M10.3d_r1", "1457-M10.3d_r2", "1457-M10.3d_r3","1457.3d_r1", "1457.3d_r2"
#)
#RNASeq.NoCellLine_reordered <- RNASeq.NoCellLine[, column_order]
#head(RNASeq.NoCellLine_reordered)
#clustering methods: "ward.D", "ward.D2", "single", "complete", "average" (= UPGMA), "mcquitty" (= WPGMA), "median" (= WPGMC) or "centroid" (= UPGMC).  pearson or spearman
datamat = RNASeq.NoCellLine[GOI, ]
write.csv(as.data.frame(datamat), file ="DEGs_heatmap_data.csv")

hr <- hclust(as.dist(1-cor(t(datamat), method="pearson")), method="complete")
hc <- hclust(as.dist(1-cor(datamat, method="spearman")), method="complete")
mycl = cutree(hr, h=max(hr$height)/1.01)
mycol = c("YELLOW", "DARKBLUE", "DARKORANGE", "DARKMAGENTA", "DARKCYAN", "DARKRED",  "MAROON", "DARKGREEN", "LIGHTBLUE", "PINK", "MAGENTA", "LIGHTCYAN","LIGHTGREEN", "BLUE", "ORANGE", "CYAN", "RED", "GREEN");
mycol = mycol[as.vector(mycl)]
png("DEGs_heatmap.png", width=1000, height=1010)
heatmap.2(as.matrix(datamat),Rowv=as.dendrogram(hr),Colv = NA, dendrogram = 'row',
            scale='row',trace='none',col=bluered(75), margins=c(4,22),
            RowSideColors = mycol, srtCol=20, keysize=0.72, cexRow = 2, cexCol = 1.4, labRow="")
dev.off()

# Generate the table for cluster members
subset_1<-names(subset(mycl, mycl == '1'))
data <- as.data.frame(datamat[rownames(datamat) %in% subset_1, ])  #212
subset_2<-names(subset(mycl, mycl == '2'))
data <- as.data.frame(datamat[rownames(datamat) %in% subset_2, ])  #185
subset_3<-names(subset(mycl, mycl == '3'))
data <- as.data.frame(datamat[rownames(datamat) %in% subset_3, ])  #173
# Take all
data <- as.data.frame(datamat)  #21859
# Initialize an empty data frame for the annotated data
annotated_data <- data.frame()
# Determine total number of genes
total_genes <- length(rownames(data))
# Loop through each gene to annotate
for (i in 1:total_genes) {
    gene <- rownames(data)[i]
    result <- getBM(attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype', 'entrezgene_id', 'chromosome_name', 'start_position', 'end_position', 'strand', 'description'),
                    filters = 'ensembl_gene_id',
                    values = gene,
                    mart = ensembl)
    # If multiple rows are returned, take the first one
    if (nrow(result) > 1) {
        result <- result[1, ]
    }
    # Check if the result is empty
    if (nrow(result) == 0) {
        result <- data.frame(ensembl_gene_id = gene,
                            external_gene_name = NA,
                            gene_biotype = NA,
                            entrezgene_id = NA,
                            chromosome_name = NA,
                            start_position = NA,
                            end_position = NA,
                            strand = NA,
                            description = NA)
    }
    # Transpose expression values
    expression_values <- t(data.frame(t(data[gene, ])))
    #colnames(expression_values) <- colnames(data)
    colnames(expression_values) <- paste(colnames(data), "normalized", sep = "_")
    # # Combine gene information and expression data
    # combined_result <- cbind(result, expression_values)
    # Fetch raw counts for the gene from count_data
    raw_counts <- count_data[gene, , drop = FALSE]
    colnames(raw_counts) <- paste(colnames(raw_counts), "raw", sep = "_")
    # Combine gene information, expression data, and raw counts
    combined_result <- cbind(result, expression_values, raw_counts)
    # Append to the final dataframe
    annotated_data <- rbind(annotated_data, combined_result)
    # Print progress every 100 genes
    if (i %% 100 == 0) {
        cat(sprintf("Processed gene %d out of %d\n", i, total_genes))
    }
}
# Save the annotated data to a new CSV file
write.csv(annotated_data, "cluster1_YELLOW.csv", row.names=FALSE)
write.csv(annotated_data, "cluster2_DARKBLUE.csv", row.names=FALSE)
write.csv(annotated_data, "cluster3_DARKORANGE.csv", row.names=FALSE)
write.csv(annotated_data, "annotated.csv", row.names=FALSE)
#~/Tools/csv2xls-0.4/csv_to_xls.py annotated.csv -d',' -o raw_and_normalized_counts.xls