add scripts

Author: QingfeiPan

.DS_Store

@@ -7,31 +7,50 @@ parent: Pipeline Setup
 
 ## Part II: Database Preparation
 
-In addition to the tools and dependencies, you will also need to prepare **databases for reference genomes** which will be used in alignment, quantification and QC assessment.
+In addition to installiing the necessary tools and dependencies, you will also need to prepare **reference genome databases** for use in alignment, quantification and QC assessment.
 
-In this pipeline, we manage the databases in **one database per reference genome** manner. Below is an overview of database preparation (use hg38 as an example):
+For this pipeline, each reference genome assembly has its own dedicated database folder. Below is an overview of the database preparation process and key inputs and outputs of each step:
 
 ![Picture](../figures/database_preparation.png)
 
 1. **Data collection**
 
-   There are **FOUR** files required in dataset preparation, and three of them can be  directly downloaded from websites:
+   There are **FOUR** files required for database preparation. Three of them can be directly downloaded from online resources:
 
-   - **<u>*annotation.gtf*</u>**: **Gene Annotation file in [GTF](https://biocorecrg.github.io/PhD_course_genomics_format_2021/gtf_format.html) (Gene Transfer Format) format**
+   - **<u>*annotation.gtf*</u>**: Gene Annotation file in [GTF](https://biocorecrg.github.io/PhD_course_genomics_format_2021/gtf_format.html) (Gene Transfer Format) format. 
 
+   - ***<u>transcriptome.fa</u>***: Transcriptome sequence file in [FASTA](https://www.ncbi.nlm.nih.gov/genbank/fastaformat/) format
 
-   - ***<u>transcriptome.fa</u>***: **Transcriptome sequence file in [FASTA](https://www.ncbi.nlm.nih.gov/genbank/fastaformat/) format**
+   - ***<u>genome.fa</u>***: Genome sequence file in [FASTA](https://www.ncbi.nlm.nih.gov/genbank/fastaformat/) format
 
-   - ***<u>genome.fa</u>***: **Genome sequence file in [FASTA](https://www.ncbi.nlm.nih.gov/genbank/fastaformat/) format**
+     ```bash
+     # locate to your conda env, change the path accordingly
+     cd /research_jude/rgs01_jude/groups/yu3grp/projects/software_JY/yu3grp/conda_env/bulkRNAseq_2025
+     
+     # create and change to the database folder
+     mkdir -p pipeline/databases/hg38/gencode.release48 # for annotation release 48 for hg38
+     cd pipeline/databases/hg38/gencode.release48
+     
+     # download the three files from https://www.gencodegenes.org/human/release_48.html
+     wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_48/gencode.v48.primary_assembly.annotation.gtf.gz
+     wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_48/gencode.v48.transcripts.fa.gz
+     wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_48/GRCh38.primary_assembly.genome.fa.gz
+     
+     # decompress them
+     gunzip *
+     
+     # rename them
+     mv gencode.v48.primary_assembly.annotation.gtf annotation.gtf
+     mv gencode.v48.transcripts.fa.gz transcripts.fa
+     mv GRCh38.primary_assembly.genome.fa.gz genome.fa
+     ```
 
-     ***<u>NOTE:</u>*** For human and mouse, we recommend [GENCODE](https://www.gencodegenes.org/) to download them, while for other species, we recommend [Ensembl](https://useast.ensembl.org/info/data/ftp/index.html). 
+     ***<u>NOTE:</u>*** For the three files listed above,  we recommend downloading them from [GENCODE](https://www.gencodegenes.org/) for human and mouse. For other species, we suggest using [Ensembl](https://useast.ensembl.org/info/data/ftp/index.html). 
 
    - ***<u>HouseKeeping gene list</u>***: the housekeeping genes defined by [this study](https://www.sciencedirect.com/science/article/pii/S0168952513000899?via%3Dihub) (N = 3804).
 
-     
+     ***<u>NOTE:</u>*** For human, the list of housekeeping genes can be downloaded [here](https://github.com/jyyulab/bulkRNAseq_quantification_pipeline/blob/main/databases/housekeepingGenes_human.txt). For other species, you can generate the housekeeping gene list by gene homology conversion using BiomaRt or other tools. Below is the codes we used to generate the housekeeping genes for mouse (avaibale [here](https://github.com/jyyulab/bulkRNAseq_quantification_pipeline/blob/main/databases/housekeepingGenes_mouse.txt)):
 
-     ***<u>NOTE:</u>*** For human, the housekeeping genes can be downloaded [here](https://github.com/jyyulab/bulkRNAseq_quantification_pipeline/blob/main/databases/housekeepingGenes_human.txt); While for other species, you can generate the housekeeping gene list by gene homology conversion using BiomaRt or other tools. Below is the codes I used to generate the housekeeping genes in mouse (avaibale [here](https://github.com/jyyulab/bulkRNAseq_quantification_pipeline/blob/main/databases/housekeepingGenes_mouse.txt)):
-     
      ``` R
      library(NetBID2)
      
@@ -47,12 +66,12 @@ In this pipeline, we manage the databases in **one database per reference genome
 
 2. **Parsing annotation file**
 
-   In this step, we will parse the gene annotation file and generate four files that required in downstream analysis:
+   In this step, we will parse the gene annotation file to generate four files required for downstream analysis:
 
-   - **`annotation.gene2transcript.txt`** & **`annotation.transcript2gene.txt`**: mappings between transcripts and genes. They are required in gene-level quantification.
-   - **`annotation.geneAnnotation.txt`** & **`annotation.transcriptAnnotation.txt`**: annotation files for genes and transcripts. They are required in the final gene expression matrix generation.
+   - **`annotation.gene2transcript.txt`** & **`annotation.transcript2gene.txt`**: These files provide mappings between transcripts and genes, which are necessary for gene-level quantification.
+   - **`annotation.geneAnnotation.txt`** & **`annotation.transcriptAnnotation.txt`**: These files contain detailed annotations for genes and transcripts, and are used in generating the final gene expression matrix.
 
-   To make the parsing analysis easiler, we created a script, **`parseAnnotation.pl`**, that you can easily generate the four files with it:
+   To simplify this process, we have provided a script, **`parseAnnotation.pl`**, which allows you to easily generate all four filest:
 
    ``` bash
    ## parse the gene anotation file

docs/1_pipeline_setup/2_database.md

@@ -0,0 +1,78 @@
+#!/usr/bin/perl
+
+use strict;
+use warnings;
+use File::Path qw(make_path);
+
+## We take the longest transcript to present each gene. The genes of which the longest transcript is less than 300nt were excluded.
+
+if (@ARGV == 3) {
+    if (-e $ARGV[0]) { unless ($ARGV[0] =~ /\.fa$/ig) { print "ERROR: The transcriptome sequence file, $ARGV[0], must be in FASTA format, and the file name much be ended with '.fa'.\n"; die; }
+    } else { print "ERROR: The transcriptome sequence file, $ARGV[0], doesn't exist. Please check and retry.\n"; die; }
+    unless (-e $ARGV[1]) { print "ERROR: The housekeeping gene file, $ARGV[1], doesn't exist. Please check and retry.\n"; die; }
+    unless (-d $ARGV[2]) { print "The output directory, $ARGV[2], doesn't exist. We are generating it...\n"; make_path($ARGV[2]); }
+} else {
+    print "ERROR: Please use THREE arguments to specify 1) transcriptome sequence file, 2) housekeeping gene file, 3) directory to save output files, respectively.\n\n\tperl prepareBins.pl xxx.transcripts.fa housekeeping_genes.xxx.txt ./bulkRNAseq/genebodyBins\n\n"; die;
+}
+
+my $input_gtf = $ARGV[0];
+my $input_hk = $ARGV[1];
+my $outout_dir = $ARGV[2];
+
+open (REF, $input_gtf) or die;
+print "We are paring $input_gtf ...\n";
+my %ref = ();
+while (<REF>) {
+    chomp;
+    next unless ($_ =~ /^>ENS/);
+    my @F = split(/\|/, $_);
+    next if ($F[6] < 300);
+    $F[0] =~ /^>(.+)/; my $isoform = $1;
+    if (exists $ref{$F[5]}) {
+        for my $x (sort keys %{$ref{$F[5]}}) {
+            if ($ref{$F[5]}{$x} > $F[6]) {
+                next;
+            }
+            else {
+                delete $ref{$F[5]}{$x};
+                $ref{$F[5]}{$isoform} = $F[6];
+            }
+        }
+    }
+    else {
+        $ref{$F[5]}{$isoform} = $F[6];
+    }
+}
+close REF;
+
+my %hk_genes = ();
+open (HK, $input_hk) or die;
+print "We are collecting housekeeping genes from $input_hk ...\n";
+while (<HK>) {
+    chomp;
+    next if ($_ =~ /^#/);
+    my @HK = split(/\t/, $_);
+    $hk_genes{$HK[0]}++;
+}
+close HK;
+
+open (OUT1, "> $outout_dir/genebodyBins_all.txt") or die;
+open (OUT2, "> $outout_dir/genebodyBins_housekeeping.txt") or die;
+print "We are writing the outputs ...\n";
+for my $x (sort keys %ref) {
+    for my $y (sort keys %{$ref{$x}}) {
+        my $length = $ref{$x}{$y};
+        my $width = $length / 100;
+
+        for my $i (0 .. 99) {
+            my $start = int($width * $i) + 1;
+            my $end = int($width * ($i + 1));
+            $end = $length if ($end > $length);
+            my $k = $i + 1;
+            print OUT1 "$y\t$start\t$end\t$x\|Bin$k\n";
+            print OUT2 "$y\t$start\t$end\t$x\|Bin$k\n" if (exists($hk_genes{$x}));
+        }
+    }
+}
+close OUT1; close OUT2;
+print "Done!\n";

scripts/.DS_Store

@@ -0,0 +1,62 @@
+#!/usr/bin/perl
+
+use strict;
+use warnings;
+
+if (@ARGV == 1) {
+    unless (-e $ARGV[0]) {
+        print "ERROR: The INPUT file doesn't exist. Please check and retry.\n"; die;
+    }
+} elsif (@ARGV < 1) {
+    print "ERROR: Please specify the gene annotation file to parse!\n\n\tperl parseAnnotation.pl xxx.annotation.gtf\n\n"; die;
+} else {
+    print "ERROR: Please use only ONE argument which specifies the gene annotation file to parse!\n\n\tperl parseAnnotation.pl xxx.annotation.gtf\n\n"; die;
+}
+
+my $f_input = $ARGV[0];
+if ($f_input =~ /\.gtf$/ig) {
+    print "We are parsing $f_input ...\n";
+
+    open (IN, $f_input) or die;
+    my %tr = (); my %gene = (); my %tr2gene = (); my %gene2tr = (); my %hash = ();
+    while (<IN>) {
+        chomp;
+        next if ($_ =~ /^#/);
+        my @F = split(/\t/, $_);
+        next unless ($F[2] eq "transcript");
+        $F[8] =~ /gene_id "(ENS.+)"; transcript_id "(ENS.+)"; gene_type "(\w+)";.* gene_name "(.+)"; transcript_type "(\w+)";.* transcript_name "(.+)"; level /;
+        $tr{$2} = "$2\t$6\t$5\t$1\t$4";
+        $gene{$1} = "$1\t$4\t$3";
+        $tr2gene{$2} = "$2\t$1"; $gene2tr{$2} = "$1\t$2"; $hash{$2} = "$2\t$6\t$1\t$4";
+    }
+    close IN;
+
+    print "We are writing the outputs ...\n";
+    (my $f_output = $f_input) =~ s/\.gtf$//g;
+    open (OUT1, "> $f_output.transcriptAnnotation.txt") or die;
+    print OUT1 "transcript_id\ttranscript_name\ttranscript_type\tgene_id\tgene_name\n";
+    for my $x (sort keys %tr) {
+        print OUT1 "$tr{$x}\n";
+    }
+    close OUT1;
+
+    open (OUT2, "> $f_output.geneAnnotation.txt") or die;
+    print OUT2 "gene_id\tgene_name\tgene_type\n";
+    for my $x (sort keys %gene) {
+        print OUT2 "$gene{$x}\n";
+    }
+    close OUT2;
+
+    open (OUT3, "> $f_output.transcript2gene.txt") or die;
+    open (OUT4, "> $f_output.gene2transcript.txt") or die;
+    for my $x (sort keys %hash) {
+        print OUT3 "$tr2gene{$x}\n";
+        print OUT4 "$gene2tr{$x}\n";
+    }
+    close OUT3; close OUT4;
+    print "Done!\n";
+}
+else {
+    print "ERROR: The INPUT file must be in GTF(Gene Transfer Format) format, and the file name much be ended with '.gtf' or '.GTF'.\n";
+}
+