cgeb2001's DokuWiki!

User Tools

Site Tools

Trace: • deconseq_for_decontaminating_reads
gene_prediction_with_braker2_pipeline

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
gene_prediction_with_braker2_pipeline [2023/02/24 18:04] – 134.190.247.138gene_prediction_with_braker2_pipeline [2025/11/18 14:23] (current) – [Genome-guided transcriptome assembly] 134.190.191.148
Line 1: Line 1:
-===== Gene prediction with Braker2 pipeline =====+====== Gene prediction with the Braker2 pipeline======
  
 GP using machine learning and extrinsic hints by **DE Salas-Leiva** (last updated Oct-21-2020)\\ GP using machine learning and extrinsic hints by **DE Salas-Leiva** (last updated Oct-21-2020)\\
Line 24: Line 24:
    mv gm_key_64 .gm_key    mv gm_key_64 .gm_key
  
-==== Repeat masking ====+===== Repeat masking ===== 
 + 
 +From the BRAKER1 paper: 
 + 
 +"Repetitive sequences create challenges for automatic gene finders both at parameter estimation step and gene prediction step. The size and quality of the training set generated by GeneMark-ET for AUGUSTUS (multi-exon genes with so called anchored introns, the introns predicted //ab initio// and also supported by RNA-Seq read mapping) is not significantly affected by TEs masking since TEs have not anchored introns. However, at the prediction step TEs can corrupt gene prediction. For this reason, soft masking of genomic sequence is recommended before execution of BRAKER1." 
 + 
 +Some repetitive elements in your genome may per-chance look like ORFs or even protein coding genes. The main purpose of masking these repeats is to prevent your gene predictor from even looking at these regions, so they will not predict any false positive genes there.
  
 Mask the repetitive regions in your assembly using the following shell script. BuildDatabase and RepeatModeler will create a species-specific library of repeats from your genome, and then RepeatMasker will use that library to mask repetitive regions in your assembly. Mask the repetitive regions in your assembly using the following shell script. BuildDatabase and RepeatModeler will create a species-specific library of repeats from your genome, and then RepeatMasker will use that library to mask repetitive regions in your assembly.
Line 71: Line 77:
 </code> </code>
  
-==== RNAseq mapping ====+===== RNAseq mapping =====
  
 +RNAseq data is direct evidence of which areas of your genome are expressed. Mapping your RNAseq data to your genome with a splice-aware mapper such as Hisat2 will yield information on the starts and stops of protein coding genes, as well as, most importantly perhaps, the start and stop coordinates of introns.
  
 On the masked assembly map the RNAseq using Hisat2, sort the output and create a depth file: On the masked assembly map the RNAseq using Hisat2, sort the output and create a depth file:
Line 86: Line 93:
 #$ -cwd #$ -cwd
 #$ -pe threaded 10 #$ -pe threaded 10
 +
 cd $PWD cd $PWD
 source activate hisat2 source activate hisat2
 +
 hisat2-build -f your_masked_genome your_genome_index hisat2-build -f your_masked_genome your_genome_index
-hisat2 -q -p  10 --phred33 --max-intronlen 1000 -k 2 -x yourgenome_index -1 READS_R1 -2 READS_R2 --rna-strandness RF -S your_masked_genome.sam+hisat2 \ 
 +    -q \ 
 +    -p 10 \ 
 +    --phred33 \ 
 +    --max-intronlen 1000 \ 
 +    -k 2 \ 
 +    -x yourgenome_index \ 
 +    -1 READS_R1 \ 
 +    -2 READS_R2 \ 
 +    --rna-strandness RF \ 
 +    -S your_masked_genome.sam 
 source deactivate source deactivate
    
 samtools view -bS -o your_genome.sambam your_masked_genome.sam samtools view -bS -o your_genome.sambam your_masked_genome.sam
-samtools sort your masked_genome.sambam -o yourmasked_genome.sambamsorted.bam+samtools sort your_masked_genome.sambam -o yourmasked_genome.sambamsorted.bam
 samtools index your_masked_genome.sambamsorted.bam samtools index your_masked_genome.sambamsorted.bam
 </code>   </code>  
          
-==== Genome-guided transcriptome assembly ====+NOTE: The ''--rna-strandness'' option here can only be used if your RNAseq experiment was done with stranded chemistry, i.e. that we know for each read whether it came from the positive or negative strand of the genome. This option is critical if we later want to split the BAM file into two BAM files, one for positive strand read mappings, and one for negative strand read mappings. The strand information in the BAM file generated by Hisat2 is captured by the ''XS'' tag. ''XS:A:+'' for positive strand mappings, and ''XS:A:-'' for negative strand mappings. If you don't invoke ''--rna-strandness'' but still perform the BAM split, for some reason only read mappings that were spliced during alignment are retained. 
 + 
 + 
 +===== Genome-guided transcriptome assembly =====
  
  
Line 107: Line 130:
 #$ -cwd #$ -cwd
 #$ -pe threaded 10 #$ -pe threaded 10
 +
 cd $PWD cd $PWD
 +
 source activate trinity-2.11-with-workaround source activate trinity-2.11-with-workaround
-Trinity --CPU 10 --max_memory 100G --genome_guided_bam yourgenome.fasta.sambamsorted.bam --genome_guided_max_intron 1000 --SS_lib_type RF+ 
 +Trinity \ 
 +    --CPU 10 \ 
 +    --max_memory 100G \ 
 +    --genome_guided_bam yourgenome.fasta.sambamsorted.bam \ 
 +    --genome_guided_max_intron 1000 \ 
 +    --SS_lib_type RF 
 conda deactivate conda deactivate
 </code> </code>
Line 122: Line 154:
  
  
-==== Braker2 ====+===== Braker2 =====
  
 +[[https://github.com/Gaius-Augustus/BRAKER|Braker]] is a fully automated pipeline in which
 +
 +  - Intron start and end coordinates (//intron hints//) are extracted from the RNAseq BAM file
 +  - These are then used along with the genome FASTA file to train GeneMarkET
 +  - The trained GeneMarkET performs an "//ab initio//" gene prediction
 +  - Those predicted gene structures for which all introns are supported by the RNAseq data (//anchored introns//) are selected to train AUGUSTUS
 +  - The trained AUGUSTUS now predicts gene structures using again the intron hints as "extrinsic evidence"
 +
 +{{::braker1_pipeline.png|}}
 +
 +The intron hints are extracted using a the ''bam2hints'' tool, with flag ''--intronsonly'', which comes with AUGUSTUS and BRAKER tools.
 +
 +If you only use RNAseq as extrinsic evidence, you essentially can only use //donor splice site// and //acceptor splice site// hints. If you also have protein homology information, you can also infer and use //start//, //stop//, //exonpart// and //exon// hints (Stanke et al 2006)
 +
 +The intron hints contain explicit location information and influence the optimal path through the GHMM machine (AUGUSTUS+ paper, Stanke et al 2006). It is important to note that since this is a probabilistic model, **hints can sometimes be ignored if the intrinsic information is strong enough!**
  
 Predict genes using Genemark-ET and Augustus through braker2: Predict genes using Genemark-ET and Augustus through braker2:
Line 152: Line 199:
  
  
-==== Predicting gene models with PASA ====+===== Predicting gene models with PASA =====
  
 PASA will use the genome-guided transcriptome assembly to estimate where gene models are located. It does this by aligning the assembled transcripts to the reference genome. PASA will use the genome-guided transcriptome assembly to estimate where gene models are located. It does this by aligning the assembled transcripts to the reference genome.
Line 206: Line 253:
 </code> </code>
  
-==== Combining gene models with EVidenceModeler (EVM) (work in progress)====+===== Combining gene models with EVidenceModeler (EVM) (work in progress) =====
  
-==NOTE: EVM requires genome sequences in Fasta format and all gene structures and alignment evidences described in GFF3 format.==+NOTE: EVM requires genome sequences in Fasta format and all gene structures and alignment evidences described in GFF3 format.==
  
-(1) Check the integrity of the gene models you wish to combine via EVM's native validator.+==== Validation ==== 
 + 
 +Check the integrity of the gene models you wish to combine via EVM's native validator.
  
 <code> <code>
Line 235: Line 284:
 </code> </code>
  
-(2) EVM requires partitioning of the genome assembly for parallelization, which greatly reduces runtime. +==== Partition ==== 
 + 
 +EVM requires partitioning of the genome assembly for parallelization, which greatly reduces runtime. 
 This step would generate directories containing fragments from each gene models for each contig, as well as a comprehensive list of these fragments and their corresponding system locations. This step would generate directories containing fragments from each gene models for each contig, as well as a comprehensive list of these fragments and their corresponding system locations.
  
Line 267: Line 318:
 </code> </code>
  
-(3) Since gene models can contain conflicting predictions, a tab-delimited weights file needs to be created to instruct EVM of which prediction to keep. +==== Assigning Weights ==== 
 + 
 +Since gene models can contain conflicting predictions, a tab-delimited weights file needs to be created to instruct EVM of which prediction to keep. 
  
 The weights file has three columns including class, type and weight. The string for class is limited to the options shown below, but other columns can have arbitray values of the same type. Weight is in ascending order. Below is an example: The weights file has three columns including class, type and weight. The string for class is limited to the options shown below, but other columns can have arbitray values of the same type. Weight is in ascending order. Below is an example:
Line 282: Line 335:
 </code> </code>
  
-(4) With all the requiste input files ready, a series of commands need to be generated for running parallelly.+==== Generating Commands ==== 
 + 
 +With all the requiste input files ready, a series of commands need to be generated for running parallelly.
  
 The file <commands.list> is generated with the below script: The file <commands.list> is generated with the below script:
Line 317: Line 372:
 </code> </code>
  
-(5) Run the series of commands in parallel with the below script:+==== Running Commands ==== 
 + 
 +Run the series of commands in parallel with the below script:
  
 <code> <code>
Line 341: Line 398:
 </code> </code>
  
-(6) Recombine the fragments for each contig:+==== Combining Fragments ==== 
 + 
 +Recombine the fragments for each contig:
  
 <code> <code>
Line 364: Line 423:
 </code> </code>
  
-(7) Convert contigs to gff3:+==== Converting To GFF3 ==== 
 + 
 +Convert contigs to gff3:
  
 <code> <code>
Line 388: Line 449:
 </code> </code>
  
-(8) Combine gff3 file from each contig to one file:+==== Combining Contigs ==== 
 + 
 +Combine gff3 file from each contig to one file:
 <code> <code>
 cat */evm.out.gff3 > EVM.all.gff3 cat */evm.out.gff3 > EVM.all.gff3
 </code> </code>
gene_prediction_with_braker2_pipeline.1677276269.txt.gz · Last modified: by 134.190.247.138