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