# EVModeler / tandy for dros spp # -------------------------------- # 0. data annotation sources: set dpid=dmoj cp -p $scd/${dpid}.fa.gz . ; gunzip ${dpid}.fa.gz echo '##gff-version 3' > ${dpid}_predict.gff # for 2nd tandy run, use limited predictor pallet: # gzcat \ $sc/caf1a/ncbi/${dp}_caf1_NCBI_GNO.gff.gz \ $sc/caf1a/pach/${dp}_CAF1.gff.gz \ $sc/caf1a/rgui/geneidv1.2/${dp}_caf1.gff3.gz \ $sc/caf1a/bren/${dp}_caf1.gff3.gz \ $sc/caf1a/batz/contrast_na/${dp}.gff.gz \ $sc/caf1c/oxfd/${dp}_caf1_oxfx.gff.gz \ $sc/caf1c/eise/glean/allsets/gff3/${dp}.gff3.gz \ $sc/caf1c/eise/glean2/gff3/${dp}.gff3.gz \ >> ${dpid}_predict.gff # recode exon >> CDS for DGIL_SNO, eise genewise, exonerate gzcat \ $sc/caf1a/dgil/${dp}_caf1_DGIL_SNO.gff.gz \ $sc/caf1a/eise/genewise/${dp}.gff3.gz \ $sc/caf1a/eise/exonerate/${dp}.gff3.gz \ | sed -e's/exon/CDS/' -e's/gene/mRNA/' >> ${dpid}_predict.gff touch ${dpid}_prothsp.gff ## recode prothsp to Target= format # Parent=CG17917_G1;tkey=CG17917-PA-HSP:46-200;tloc=46-200;align=155 gzcat \ $sc/caf1a/dgil/${dp}prot9-hsp.gff.gz \ | perl -pe's/HSP/protein_match/; s/tkey=([\w\-:]+)-HSP:(\d+)-(\d+)/Target=$1 $2 $3/;' \ >> ${dpid}_prothsp.gff # maybe later: # $sc/caf1a/dgil/${dp}prot9-hsp.gff.gz \ << HSP # $sc/caf1c/dgil/${dp}_caf1_DGIL_TEX.gff.gz \ << match # # $sc/caf1a/eise/genemapper/${dp}.gff3.gz \ # $sc/caf1a/oxfd/${dp}_caf1.gff.gz \ # $sc/caf1a/dgil/${dp}_caf1_DGIL.gff.gz \ ------------ # 1. partition # ** REVISE THIS, for tandy use, split at 2 MB segs with 50Kb overlap, min 1 MB segs $EVM_HOME/EvmUtils/partition_EVM_inputs.pl \ --genome ${dpid}.fa \ --gene_predictions ${dpid}_predict.gff \ --protein_alignments ${dpid}_prothsp.gff \ --segmentSize 500000 --overlapSize 10000 \ --minSegmentSize 1000000 \ --partition_listing ${dpid}_partitions.list \ #>> fixme: --segmentSize 500000 --overlapSize 10000 \ # --transcript_alignments ${dpid}_estmatch.gff \ # --pasaTerminalExons ${dpid}_terminal_exons.gff \ #---- set dp=dmoj set species=mojavensis set dpid=${dp}_caf060210 set bmid=61 set scd=$sc/${dp}3 set dp3=${dp}3 #------------- prelim. tandy 5n test on small scaffold # run2 limited predict set: leave out GI_PACH_GMP_, GI_EISE_CEX_, GI_BATZ_CNA_, TRdmoj_?? == oxfd, # keep GI_BREN_NSC_, GI_DGIL_SNO_, GI_NCBI_GNO_, GLEAN_, dmoj_GLEANR_?? #------ Daphnia 4.e ------------- # 2007jun21 # 0. data annotation sources: set dpid=dpulex1 # genome.fa: hard mask the softmask in this genome cp -p $sc/dpulex4/dappu1jgi/dpulex_jgi060905.repeatmasked.fa.gz ${dpid}.fa.gz gunzip ${dpid}.fa.gz perl -pi -e'unless(/^>/){s/[a-z]/n/g;}' ${dpid}.fa ## gffs: predict, prothsp, prot4gw [1-9], estmatch, goldgene ? from pasa echo '##gff-version 3' > ${dpid}_predict.gff gzcat \ $sc/dpulex4/dpulex_jgi060905_Gnomon.gff.gz \ >> ${dpid}_predict.gff # remove DGIL_SNO from JGI predicts gzcat $sc/dpulex4/dpulex_jgi060905_JGI_FM5.gff.gz \ | perl -ne 'if(/SNAP/) { ($sd)= m/(Dappu1_FM5_\d+)/; $snap{$sd}++; next; }\ elsif(/exon|CDS/){ ($d)=m/Parent=(\w+)/; next if($snap{$d}); } print;'\ > ${dpid}_jginospan.gff cat ${dpid}_jginospan.gff >> ${dpid}_predict.gff # recode exon >> CDS for DGIL_SNO, eise genewise, exonerate gzcat \ $sc/dpulex4/dpulex_jgi060905_DGIL_SNO.gff.gz \ | sed -e's/exon/CDS/' -e's/gene/mRNA/' \ > ${dpid}_DGIL_SNO.gff cat ${dpid}_DGIL_SNO.gff >> ${dpid}_predict.gff ## recode prothsp to Target= format touch ${dpid}_prothsp.gff gzcat \ $sc/dpulex4/dpulex-prot4-hsp.gff.gz \ | perl -pe's/HSP/protein_match/; s/tkey=([\w\-:]+)-HSP:(\d+)-(\d+)/Target=$1 $2 $3/;' \ >> ${dpid}_prothsp.gff ## pasa ESTs, partial GeneWise (scaff 1..9?) cp -p $sc/dpulex4/pasa_daphc.pasa_assemblies.gff.gz ${dpid}_estmatch.gff.gz gunzip ${dpid}_estmatch.gff.gz # cp -p ../daphd/dpulex1_prot4gw.gff.gz . ; gunzip dpulex1_prot4gw.gff.gz # goldgene = pasa_daphc.genemodel_updates.all.gff.gz ? << full models only $em/PasaUtils/retrieve_terminal_CDS_exons.pl \ $pa/daphc/trainingSetCandidates.fasta $pa/daphc/trainingSetCandidates.gff \ > ! ${dpid}_terminal_exons.gff # 1. partition # ** remove bacteria scaffolds, ?? drop small segs < 10kb? $EVM_HOME/EvmUtils/partition_EVM_inputs.pl \ --genome ${dpid}.fa \ --gene_predictions ${dpid}_predict.gff \ --protein_alignments ${dpid}_prothsp.gff \ --transcript_alignments ${dpid}_estmatch.gff \ --pasaTerminalExons ${dpid}_terminal_exons.gff \ --segmentSize 1000000 --overlapSize 50000 \ --minSegmentSize 10000 \ --skipSegmentList dpulex_jgi060905_remove_bact_scaffs1.list \ --partition_listing ${dpid}_partitions.list \ > & log.parts & melon.% grep -c 'skipped' log.parts 7874 melon.% ls -d1 scaffold* | wc 1206 melon.% grep -c ' has ' log.parts 1206 #------- Celegans test ----------- # 0. data annotation sources: set dpid=cele1 Use this genome data, 7 chromosomes (1 is MtDNA) $sc/cele1/celegans-dna-WS167.fa.gz microbe% zgrep mRNA $sc/cele1/celegans-slim-WS167.gff.gz | perl -ne'($r,$s,$t,@x)=split; print "$s\n" if($t eq "mRNA");' | sort | uniq -c 27049 Coding_transcript << these mRNA microbe% zgrep Coding_transcript $sc/cele1/celegans-slim-WS167.gff.gz | perl -ne'($r,$s,$t,@x)=split; print "$t\n";' | sort | uniq -c 171044 CDS << tandy data set, why more than exons? 136862 exon 17515 five_prime_UTR 20089 gene 27049 mRNA 14810 three_prime_UTR #................. set dpid=cele1 cp -p $sc/cele1/celegans-dna-WS167.fa.gz ${dpid}.fa.gz ; gunzip ${dpid}.fa.gz ## worm gff needs (a) recoding or (b) adjust tandemgenes.pl ## .. remove ID= from CDS lines, using only Parent= ## .. remove leading 'Transcript:' etc tags from Parent/ID values (why there??) ## all wormbase IDs seem to have a type prefix (Gene:,Transcript:,CDS:) ## ** ALSO bug for tandemgenes, worm IDs have digit suffixes, confusing tandy's exon digit additions: ## recode . to _ for all ids : turning all attrib '.' to _ is ugly but should work echo '##gff-version 3' > ${dpid}_predict.gff gzgrep Coding_transcript $sc/cele1/celegans-slim-WS167.gff.gz | perl -pe \ '@v=split"\t"; s/ID=[^;]+;// if($v[2] =~ /(CDS|exon|UTR)/); s/(ID|Parent)=\w+:/$1=/g; \ @v=split"\t"; $v[8] =~ s/\./_/g; $_= join("\t",@v); ' \ >> ${dpid}_predict.gff touch ${dpid}_prothsp.gff gzcat \ $sc/cele1/celeprot9-hsp.gff.gz \ | perl -pe's/HSP/protein_match/; s/tkey=([\w\-:]+)-HSP:(\d+)-(\d+)/Target=$1 $2 $3/;' \ >> ${dpid}_prothsp.gff $EVM_HOME/EvmUtils/partition_EVM_inputs.pl \ --genome ${dpid}.fa \ --gene_predictions ${dpid}_predict.gff \ --protein_alignments ${dpid}_prothsp.gff \ --segmentSize 1000000 --overlapSize 50000 \ --partition_listing ${dpid}_partitions.list \ >& parts.log & # --minSegmentSize 1000000 \ NOTE: CDS/mRNA ids are not predictor-prefixed, so stats like this don't work: -- need some flag in tandy to say if input IDs have method prefixes... cat scaffold_*/*_exons_tandy.gff | grep bestids= | perl -n $td/genebest.perl #................. # Dros. mel. tandy 1; using only ref genes/cds # FIXME: problem with ncRNA-exons mixed in with CDS-exons # FIXME 2: really need to mask out Transposons, causing lots of extra exon matches # 0. data annotation sources: set dpid=dmel4 melon.% ls $sc/dmel3 dmel-chromosomes.gff dmel_r420shred.fa.gz@ gff@ perchr@ dmel-markers.gff.gz@ dmel_r430.fa.gz@ noexonchr/ pgsql@ dmel_r4.3_20060217@ dpulex-dmel3noexonte.gff noexontechr/ cp -p $sc/dmel3/dmel_r430.fa.gz ${dpid}.fa.gz ; gunzip ${dpid}.fa.gz echo '##gff-version 3' > ${dpid}_predict.gff ## ** FIXME: this includes non-coding RNA exons, but excludes the ncRNA/tRNA/snRNA/.. features ## we don't want any ncrna exons, see dmel.ncrna.idlist (FBtr ids) to exclude; ## the ncrna are over-abundant in genome repeats compared to CDS-exons (see esp. Dmel.chr2R) gzcat $sc/dmel3/gff/dmel-*-r4.3.0.gff.gz | grep FlyBase | perl -ne \ '($r,$s,$t)=split"\t"; print if($t =~ /gene|mRNA|exon|CDS|protein|chromosome_arm/ && $s eq "FlyBase"); ' \ >> ${dpid}_predict.gff touch ${dpid}_prothsp.gff gzcat \ $sc/caf1a/dgil/dmelprot9-hsp.gff.gz \ | perl -pe's/HSP/protein_match/; s/tkey=([\w\-:]+)-HSP:(\d+)-(\d+)/Target=$1 $2 $3/;' \ >> ${dpid}_prothsp.gff $EVM_HOME/EvmUtils/partition_EVM_inputs.pl \ --genome ${dpid}.fa \ --gene_predictions ${dpid}_predict.gff \ --protein_alignments ${dpid}_prothsp.gff \ --segmentSize 1000000 --overlapSize 50000 \ --partition_listing ${dpid}_partitions.list \ >& parts.log & #................. # Dros. mel. tandy 2; including CAF1 predictions and ref genes; excluding ncRNA # FIXME 2: really need to mask out Transposons, causing lots of extra exon matches # use $sc/dmel3/dmel4.fa.temasked.gz :: SOFTmasked; use blat -mask=lower or hardmask set dpid=dmel4c ln -s ../dmel4/dmel4.fa $dpid.fa ## same as ../dmel4/*.fa # cp -p $sc/dmel3/dmel_r430.fa.gz ${dpid}.fa.gz ; gunzip ${dpid}.fa.gz echo '##gff-version 3' > ${dpid}_predict.gff gzcat $sc/dmel3/gff/dmel-*-r4.3.0.gff.gz | grep FlyBase | grep -v FlyBase.match | \ perl -ne '($r,$s,$t)=split"\t"; ($id)=m/ID=(\w+)/; \ print "$id\n" if($t =~ /RNA/ && $t ne "mRNA" && $s eq "FlyBase"); ' \ | sort | uniq > dmel.ncrna.idlist gzcat $sc/dmel3/gff/dmel-*-r4.3.0.gff.gz | grep FlyBase | grep -v FlyBase.match | \ perl -ne \ 'BEGIN{ open(F,"dmel.ncrna.idlist"); while(){ chomp; $ncrna{$_}++;} close(F);} \ ($r,$s,$t)=split"\t"; m/(ID|Parent)=(\w+)/; $id=$2; $skipnc++ and next if($ncrna{$id}); \ END{warn"skipped $skipnc ncRNA features\n";}\ print if($t =~ /gene|mRNA|exon|CDS|protein|chromosome_arm/ && $s eq "FlyBase"); ' \ >> ${dpid}_predict.gff ^^ *** this made mistake and left in dmel.ncrna.idlist items ; bad read list*** # ? collect these: transposable_element to check if predictor exon matches # are hitting transposon sites? gzcat $sc/dmel3/gff/dmel-*-r4.3.0.gff.gz | grep FlyBase | grep FlyBase.transpos \ > ${dpid}_transposon.gff ## problem: CDS vs exon; above dmel ref has only exon + protein range to subset ## need to exclude exon** from predictors, or exclude CDS # note: batz/contrast has start/stop_codon + CDS; GNO has gene/mRNA/exon,CDS; # rgui has gene/mRNA/CDS ; oxfd has gene/mRNA/CDS gzcat \ $sc/caf1a/ncbi/${dp}_caf1_NCBI_GNO.gff.gz \ $sc/caf1a/rgui/geneidv1.2/${dp}_caf1.gff3.gz \ $sc/caf1a/batz/contrast/${dp}.gff.gz \ $sc/caf1c/oxfd/${dp}_caf1_oxfx.gff.gz \ | grep -v 'NCBI_GNO.exon' >> ${dpid}_predict.gff # recode exon >> CDS for DGIL_SNO gzcat \ $sc/caf1a/dgil/${dp}_caf1_DGIL_SNO.gff.gz \ | sed -e's/exon/CDS/' -e's/gene/mRNA/' >> ${dpid}_predict.gff ## feature counts # 54058 CDS BATZ_CON # 79074 CDS DGIL_SNO # 82949 CDS NCBI_GNO # 84543 CDS OXFD_GPX # 53111 CDS RGUI_GID # 8 chromosome_arm FlyBase # 63346 exon FlyBase # 14486 gene FlyBase # 15538 gene NCBI_GNO # 13572 gene OXFD_GPX # 12671 gene RGUI_GID # 23937 mRNA DGIL_SNO # 19389 mRNA FlyBase # 20420 mRNA NCBI_GNO # 19115 mRNA OXFD_GPX # 12671 mRNA RGUI_GID # 19389 protein FlyBase # 1372 protein_binding_site FlyBase # 51 pseudogene FlyBase # 14090 start_codon BATZ_CON # 14103 stop_codon BATZ_CON ## feature counts version with ref ncrna mistake: # 64150 exon FlyBase ** has 1000 ncrna exons # 1. partition # ** REVISE THIS, for tandy use, split at 2 MB segs with 50Kb overlap, min 1 MB segs $EVM_HOME/EvmUtils/partition_EVM_inputs.pl \ --genome ${dpid}.fa \ --gene_predictions ${dpid}_predict.gff \ --segmentSize 1000000 --overlapSize 50000 \ --minSegmentSize 1000000 \ --partition_listing ${dpid}_partitions.list \ > & parts.log & #later: --protein_alignments ${dpid}_prothsp.gff \