bash ./module-check.sh

bash ./ref-download.sh -o mouse -r M22 -t 4

bash ./orf-prediction.sh -o \"Mus musculus\" -t 8

bash ./riboseq-process.sh -f ./out/data/ribo-seq/ribo.fastq.gz -a AAAAAAAAAAA -t 4

bash ./rnaseq-process.sh -f ./out/data/rna-seq/rna.fastq.gz -t 4

bash ./orf-calling.sh -o mouse -x 10090 -m 32 -n 28 -t 8

# Tools to install on CentOS before miniconda
yum -y install gcc tar bzip2 git which

curl -fsSL https://repo.anaconda.com/miniconda/Miniconda2-latest-Linux-x86_64.sh -o miniconda2.sh

# assume miniconda is installed in the home directory
bash miniconda2.sh -b -p ~/miniconda2

export PATH=~/miniconda2/bin:$PATH
export PYTHONPATH=~/miniconda2/lib/python2.7/site-packages

conda install -y -c conda-forge wget 
conda install -y -c conda-forge parallel
conda install -y -c bioconda samtools
conda install -y -c bioconda htslib 
conda install -y -c bioconda bedtools 
conda install -y -c bioconda bedops 
conda install -y -c bioconda bowtie 
conda install -y -c bioconda fastqc
conda install -y -c bioconda cutadapt 
conda install -y -c bioconda trim-galore 
conda install -y -c bioconda star 
conda install -y -c bioconda stringtie 
conda install -y -c bioconda sra-tools
conda install -y -c bioconda emboss 
conda install -y -c bioconda plastid
conda install -y -c bioconda bioconductor-rhtslib
Rscript -e 'install.packages("BiocManager", repos="https://cran.us.r-project.org"); BiocManager::install(c("riboSeqR", "GenomicFeatures", "rtracklayer"))'

git clone https://github.com/boboppie/ORFLine.git
cd ORFLine
chmod +x *.sh
  
bash ./main.sh

singularity pull shub:https://boboppie/ORFLine-singularity # the default name of the image is ORFLine-singularity_latest.sif

singularity run ORFLine-singularity_latest.sif

# Assume the source code is in ~/code/github/ORFLine/
# Mouse GENCODE M13 reference files are in ~/ref/GENCODE/mouse/M13/
# Transcript sequence file downloaded from GENCODE is named gencode.vM13.transcripts.fa.gz, and placed in ~/ref/GENCODE/mouse/M13/fasta/transcriptome/
# The newly generated transcript sequence file name is transcripts_biotype_filtered.fa

cd ~/ref/GENCODE/mouse/M13/fasta/transcriptome/

~/code/github/ORFLine/script/fastagrep.sh -v 'IG_*|TR_*|miRNA|misc_RNA|Mt_*|rRNA|scaRNA|scRNA|snoRNA|snRNA|sRNA|ribozyme|nonsense_mediated_decay|non_stop_decay' <(gzip -dc <gencode_transcripts.fa.gz>) > <transcriptome.fa>

# Arguments:
# START_CODON - start codon, can be “ATG”, “TTG”, “CTG”, “GTG”  
# TRANSCRIPTOME_FASTA - transcriptome file location, assume is it in ~/ref/GENCODE/mouse/M13/fasta/transcriptome
# GTF - gene annotation file, assume it is in ~ref/GENCODE/mouse/M13/annotation/CHR/comprehensive/, and named gencode.vM13.annotation.gtf
# ORGANISM - scientific name, e.g. Mus musculus
# NCORE - number of computing cores if you want the process to be multi-threading, default 1
#
# for example, they can be assigned as: 
# START_CODON=ATG
# TRANSCRIPTOME_FASTA=~/ref/GENCODE/mouse/M13/fasta/transcriptome/transcripts_biotype_filtered.fa
# GTF=~ref/GENCODE/mouse/M13/annotation/CHR/comprehensive/gencode.vM13.annotation.gtf
# ORGANISM="Mus musculus"
# NCORE=8 
#
# Output:
# The output file is in BED12 format (https://genome.ucsc.edu/FAQ/FAQformat#format1) and named orf_$START_CODON.bed 

Rscript ~/code/github/ORFLine/script/ORFPredict.R <start_codon> <transcriptome.fa> <annotation.gtf> <organism_scientific_name> <threads> 

cut -f11 <orfs_ATG.bed> | sed -e 's/,/\t/g' | awk '{for(i=t=0;i<NF;) t+=$++i; $0=t}1' > <orfs_ATG.width>

paste -d'\t' <orfs_ATG.bed> <orfs_ATG.width> | awk '$13 <= 303' | cut -f1-12 | sort -k4 > <orfs_ATG.smORFs_100.bed>

fastqc -t <threads> -o <fastqc_dir> <sample.fastq.gz>

# User can give a customized adapter by flag -a, or trim_galore will automatically detect the adapter if it's standard
# To filter for read length, user can use flags such as --length 25 --max_length 35 for min and max length

trim_galore -q 33 --fastqc --trim-n -e 0.1 --stringency 3 <sample.fastq.gz> -o <trim_galore_dir>

# All contaminant sequeneces can be merged in a file, assume named contaimination.fa

bowtie-build <contaimination.fa> <contaimination_bowtie_index_dir> 

bowtie -a --best --strata -S --seed 23 -p <threads> --chunkmbs 256 --norc --maqerr=60 --un <trimmed_unfiltered.fq> <contaimination_bowtie_index_dir> <(gzip -dc <trimmed.fq.gz>) <trimmed_filtered.sam>

gzip <trimmed_unfiltered.fq>

# sjdbOverhang is the longest read length - 1

STAR --runThreadN <threads> --runMode genomeGenerate --genomeDir <star_ribo_dir> --genomeFastaFiles <ref_genome.fa> --sjdbGTFfile <annotation.gtf> --sjdbOverhang <sjdbOverhang>

SAM_ATTR="NH HI NM MD AS"
ALIGNINTRON_MIN=20
ALIGNINTRON_MAX=10000
MISMATCH_MAX=1
MISMATCH_NOVERL_MAX=0.04
FILTER_TYPE=BySJout

STAR --runThreadN <threads> --genomeDir <star_ribo_dir> \
     --readFilesIn <trimmed_unfiltered.fq.gz> --readFilesCommand zcat \
     --outReadsUnmapped Fastx --outFileNamePrefix <prefix> \
     --alignIntronMin <alignIntronMin> --alignIntronMax <alignIntronMax> --alignEndsType EndToEnd \
     --outFilterMismatchNmax <mismatch_max> --outFilterMismatchNoverLmax <mismatch_noverl_max>\
     --outFilterType <filter_type> --outFilterIntronMotifs RemoveNoncanonicalUnannotated \
     --outSAMattributes <sam_attr> --outSAMtype BAM SortedByCoordinate --outBAMsortingThreadN <threads>
     
# Only uniquely mapped reads are kept for ORF calling

samtools view -b -q 255 <Aligned.sortedByCoord.out.bam> > <sample_q255.bam>
samtools index <sample_q255.bam>

# We keep protein coding leve 1 and 2, not seleno proteins

cat <gencode_annotation.gtf> | awk '{if($18=="\"protein_coding\";" && $0~"level (1|2);" && $0!~"tag \"seleno\";" && $0!~"cds_end_NF" && $0!~"cds_start_NF"){print $0}}' | sort -k1,1 -k4,4n | bgzip > <protein_coding.gtf.gz>

tabix -p gff <protein_coding.gtf.gz>

metagene generate <metagene_analysis_dir> \
                  --landmark cds_start \
                  --annotation_files <protein_coding.gtf.gz>

psite <metagene_analysis_dir>/protein_coding_rois.txt <psite_dir> \
      --min_length 25 \
      --max_length 35 \
      --require_upstream \
      --count_files <sample_q255.bam>

phase_by_size <metagene_analysis_dir>/protein_coding_rois.txt <phasing_dir> \
              --count_files <sample_q255.bam> \
              --fiveprime_variable \
              --offset <psite_dir>/p_offsets.txt \
              --codon_buffer 5

nice -5 STAR --runThreadN <threads> --genomeDir <star_rna_100_pe_dir> \
                       --readFilesIn <R1_trimmed.fq.gz> <R2_trimmed.fq.gz> --readFilesCommand zcat \
                       --alignEndsType EndToEnd \
                       --outFilterMismatchNmax <mismatch_max> \
                       --outReadsUnmapped Fastx --outFileNamePrefix <prefix> \
                       --outSAMattributes <sam_attr> --outSAMunmapped Within --outSAMtype BAM SortedByCoordinate \
                       --outBAMsortingThreadN <threads>

stringtie <sample_q255.bam>  -p <threads> -G <protein_coding.gtf> -eB -o expressed.gtf -A gene_abund.tab -C cov_refs.gtf

samtools view -h <ribo.bam> | awk -v minl="<min_len>" -v maxl="<max_len>" 'length($10) >= minl && length($10) <= maxl || $1 ~ /^@/' | samtools view -bS - > filtered.bam

samtools index filtered.bam

parallel -j<threads> "bedtools coverage -counts -split -a <orfs_{}.smORFs.bed> -b <filtered.bam> > <orfs_{}.smORFs.bedtools.counts>" ::: ATG CTG TTG GTG

parallel "awk '\$13 > 0' <orfs_{}.smORFs.bedtools.counts> > <orfs_{}.smORFs.bedtools.greaterThanZeroReads.counts>" ::: ATG CTG TTG GTG
parallel "cut -f1-12 <orfs_{}.smORFs.bedtools.greaterThanZeroReads.counts> > <orfs_{}.smORFs.bedtools.greaterThanZeroReads.bed>" ::: ATG CTG TTG GTG

parallel "python script/get_count_vectors.py --annotation_files <orfs_{}.smORFs.bedtools.greaterThanZeroReads.bed> --annotation_format BED --count_files <filtered.bam> --min_length <min_len> --max_length <max_len> --fiveprime_variable --offset <psite_offset_file> --out_prefix {} ." ::: ATG CTG TTG GTG

parallel "LC_ALL=C fgrep -f <stringytie_active_transcript_Ids.txt> <plastid_ORFs_{}_translated.txt> > <ORFs_{}_translated_expressed.txt>" ::: ATG CTG TTG GTG

parallel "Rscript <ORFScore.R> <ORFs_{}_translated_expressed.txt> {}" ::: ATG CTG TTG GTG

cat <ORFScore_*TG.tsv> | sort -k2 -g | awk '{if($2!="NA") print}' > <ORFScore_all_noNA.tsv>
cat <ORFScore_all_noNA.tsv> | awk '{if($2 >0 && $6 >= 0.1) print}' > <ORFScore_all_filteredByCov.tsv>

Rscript ~/code/github/ORFLine/script/RegionFilter.R . <cds.gtf> <filtered.bam> <threads>

Rscript script/RRS.R . <rna-seq_bam_dir>

tail -n +2 <RRS_out.tsv> | awk '$2 > 9 && $3 < 30 && $4 > 5' > <RRS_out_filtered.tsv>

Rscript script/NestedRegionFilter.R . <threads>

Rscript script/ORFScore_padj.R <ORFScore_merged_PT_filtered.tsv> <all_withQval>

awk '($4+0) < 0.01' <ORFScore_all_withQval.tsv> > <ORFScore_fdr_filtered.tsv>

# Length filter >= 11aa including stop codon
cat <ORFScore_fdr_filtered.tsv> | cut -f1 | cut -f4,5 -d: | tr ":" "\t" | awk '{print ($2-$1+1)/3}' > <ORFScore_fdr_filtered_len.txt>
paste <ORFScore_fdr_filtered.tsv> <ORFScore_fdr_filtered_len.txt> | awk '$17 >= 11' | cut -f1-16 > <ORFScore_length_filtered.tsv>

# codon_1_2 mean should be greater than 0
awk '$10 > 0' ORFScore_length_filtered.tsv >ORFScore_codon_filtered.tsv