Unveiling Difference in Promoter regions of MeDIP data 2

1. 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
   

2. 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 "protein_coding" tss_.bed > tss.bed #21859 records

3. 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.

4. 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
   

5. 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
   

6. 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
   

like unlike

点赞本文的读者

还没有人对此文章表态

本文有评论

没有评论

看文章，发评论，不要沉默