Asm-RNA annotation and selection steps 1 - add protein, CDS (best orf), check/correct with valid introns 2 - select best transcripts/locus 1a - add CDS, protein to asmrna transcript.gff (or check them if from gmap) ** Corrected incut-pro method : replace makeincutpro.sh ; add utrx= attrib setenv PERL5LIB ${PERL5LIB}:$workd/PerlLib/ $evigene/scripts/genefindcds.pl -dna $workd/genome/nasvit1asm.fa -intron $workd/intron/intron_good.gff.gz \ -genes $geneset.gff.gz > $geneset.pinfix.gff 1b - score for valid introns, homo proteins: requires evidence of $workd/prot/protein.gff $workd/intron/intron.gff scripts/inmatescore.sh > asmrna.{introntab,inprot} env DOAN=1 workd=$workd genes=$gset.pinfix.gff $workd/scripts/inmatescore.sh 1c - score to an1.gff (now part of inmatescore.sh) cat $gset.introntab | cut -f1,3 | sed 's,/.*,,; s/^/intr /' > $gset.intr.tmp cat $gset.inpro | cut -f1,2 | sed 's,/.*,,; s/^/pro /;' > $gset.inpro.tmp cat $gset.{intr,inpro}.tmp $gset.incutcd.gff | perl -ne \ 'if(/^(intr|pro)\s/){ ($k,$d,$v)=split; $dv{$d}{$k}=$v; } \ else{ if(/\tmRNA/){ ($d)=m/ID=([^;\s]+)/; foreach $k (qw(intr pro)) \ { if( $v=$dv{$d}{$k}) { s/$/;$k=$v/; } } } print; }' \ > $gset.an1.gff 2a - mark alt-tr (and best-tr) per locus per asmrna set, method needs work, 2a.a - overgenedup -in genes -over genes | .. > tabalt : table of gene,alt,ident : not best result? 2a.b - overbestgene2 -alttr ; some but not best method 2b - best tr/locus of each set (preserving alttr) cat $gset.an1.gff | \ $workd/scripts/overbestgene2.perl -in stdin \ -alttr -typeover CDS -mrna mRNA -exon 'CDS,exon' \ -scoretype='intr:90,pro:8,CDS:3,UTR:1' \ -drops='intr:1,pro:30,CDS:60,UTR:180' \ -OVEREXON2 -genescore -trivial 10 -summarize -skip \ > $gset.best1.gff 2c - best among sets : add option overbestgene2 : -resource/reid to change source and prefix ID/Parent tags -resource=nv1bof4 -preid=b4 cat nvit1_rnaseq*.best1.gff | \ $workd/scripts/overbestgene2.perl -in stdin \ -typeover CDS -mrna mRNA -exon 'CDS,exon' \ -scoretype='intr:90,pro:6,CDS:3,UTR:1' \ -drops='intr:1,pro:30,CDS:60,UTR:180' \ -OVEREXON2 -genescore -trivial 10 -summarize -skip \ > nvit1_rnaseq.bestof4.gff # for map views update source/ID perl -pi -e'if(/^\w/){ s/\t\w+/\tnv1bof4/; s/(ID|Parent)=/$1=b4/; s/Parent=/gene=/ if(/\tmRNA/);}' \ nvit1_rnaseq.bestof4.gff perl -pi -e'if(/^\w/){ s/\t(\w)/\tb1$1/; s/(ID|Parent)=/$1=b1/; s/Parent=/gene=/ if(/\tmRNA/);}' \ nvit1_rnaseq.*.best1.gff #....... quality checks ......... Gene assemblies coding / transcript size summary ------------------------------------------------- est-cgbasm: -- nonoverlapped (no alttr) ALL n=58207; aveCDS= 636; aveEXON=1035; %C/X=61; sumC=37028377; sumX=60277082 TOP n=25000; aveCDS=1109; aveEXON=1456; %C/X=76; sumC=27746172; sumX=36422416 -- all asm (+alttr) ALL n=96702; aveCDS=638; aveEXON=1022; %C/X=62; sumC=61709513; sumX=98865949 TOP n=25000; aveCDS=1284; aveEXON=1646; %C/X=78; sumC=32124402; sumX=41169840 rna-cufasm: -- nonoverlapped (no alttr) ALL n=18447; aveCDS=652; aveEXON=971; %C/X=67; sumC=12035624; sumX=17924030 TOP n=18447; aveCDS=652; aveEXON=971; %C/X=67; sumC=12035624; sumX=17924030 -- all asm (+alttr) ALL n=21513; aveCDS=685; aveEXON=1011; %C/X=67; sumC=14738680; sumX=21763635 TOP n=21513; aveCDS=685; aveEXON=1011; %C/X=67; sumC=14738680; sumX=21763635 estrna-velasm: -- nonoverlapped (no alttr) ALL n=152279; aveCDS= 292; aveEXON= 591; %C/X=49; sumC=44513794; sumX=90120950 TOP n=25000; aveCDS=1131; aveEXON=1859; %C/X=60; sumC=28279040; sumX=46486771 -- all asm (+alttr) ALL n=205854; aveCDS=363; aveEXON=726; %C/X=50; sumC=74830793; sumX=149653086 TOP n=25000; aveCDS=1628; aveEXON=2482; %C/X=65; sumC=40708883; sumX=62068326 cacao9consensus1 (no alttr, excluding TE genes): ALL n=28490; aveCDS=1283; aveEXON=2905; %C/X=44; sumC=36575569; sumX=82791713 TOP n=25000; aveCDS=1430; aveEXON=3187; %C/X=44; sumC=35761102; sumX=79684628 --------- Same CDS span as cacao9consensus1 est-cgbasm : 5065 << may be highest because this est set informed cons1 gene set rna-cufasm : 1786 estrna-velasm : 2793 ------------------------------------------------- cat xxx-mars11.pinfix.gff | grep -v '^#' | ggrep -F -f xxx-mars11.bestids - | perl -ne \ '($d)=m/Parent=([^;\s]+)/; \ if($ld ne $d and $ld) { $n++; $cs{$ld}= $c; $xs{$ld}=$x; $cs+=$c; $xs+=$x; $x=$c=$nx=0; } \ ($b,$e)=(split"\t")[3,4]; $w=1+$e-$b; \ if(/\tCDS/){ $c+=$w; } elsif(/\texon/) { $x+=$w; } $ld=$d;\ END{ $ca=int($cs/$n); $xa=int($xs/$n); $pc=int(100*$cs/$xs); \ print "ALL n=$n; aveCDS=$ca; aveEXON=$xa; %C/X=$pc; sumC=$cs; sumX=$xs\n"; \ @d= sort{$cs{$b}<=>$cs{$a}} keys %cs; $dm=(@d>25000)?25000:@d; @d=splice(@d,0,$dm); \ $cs=$xs=$n=0; foreach $d (@d) { $n++; $cs+=$cs{$d}; $xs+=$xs{$d}; } \ $ca=int($cs/$n); $xa=int($xs/$n); $pc=int(100*$cs/$xs); \ print "TOP n=$n; aveCDS=$ca; aveEXON=$xa; %C/X=$pc; sumC=$cs; sumX=$xs\n"; } ' flamingo2.% echo $gset ; grep mRNA $gset | grep -v '^#' | grep -v 'alttr=' | cut -f9 | perl -ne 'for $i (0..3) { ($sv,$x v)= (m/$kv[$i]=(\d+).(\d*)/)?($1,$2):(0,0); $n[$i]++ if($sv>0); $s[$i]+=$sv; $t[$i]+=$xv; } $n++; BEGIN{@kv=qw(intronfix aalen cxlen utrx); } \ END{ for $i (0..3){ $s=$s[$i]; $t=$t[$i]; $ni=$n[$i]; $at=int(100*$t/$n)/100; $a=int(100*$s/$n)/100; print join("\t",$kv [$i],$a,$s,$at,$t,$ni,$n),"\n"; }}' rscuff.cacao11rs54.pinfix.gff ; ntr=21513 intronfix 0.02 464 0 0 386 21513 aalen 227.51 4894480 66.84 1438027 21513 21513 cxlen 682.53 14683440 1012.78 21788093 21513 21513 utrx 0.07 1573 0.07 1572 417 21513 vs cacao9_consensus1_note: aalen 415.37 11837038 60.63 1728005 28496 28497 rscuff.cgb_ca9.pinfix.gff : ntr=33258; more genes, but shorter/lower qual, for 1 of 54 bams intronfix 0.02 802 0 0 690 33258 aalen 169.52 5638113 63.61 2115740 33258 33258 cxlen 508.57 16914339 802.46 26688297 33258 33258 utrx 0.05 1823 0.04 1481 486 33258 cat rscuff.cacao11rs54.pinfix.gff | grep -v '^#' | egrep ' (exon|CDS)' | perl -ne \ '($d)=m/Parent=(\w+)/; next if($d=~m/p\d+$/); if($ld ne $d and $ld) { $n++; $cs+=$c; $xs+=$x; $x=$c=$nx=0;} \ ($b,$e)=(split)[3,4]; $w=1+$e-$b; if(/CDS/){ $c+=$w; } elsif(/exon/) { $x+=$w; } $ld=$d; \ END{ $ca=int($cs/$n); $xa=int($xs/$n); $pc=int(100*$cs/$xs); \ print "n=$n; aCDS=$ca; aEX=$xa; %C/X=$pc; sC=$cs; sX=$xs\n"; } ' # rnaseq-cufflinks of 54 groups: cacao11rs54.pinfix: n=21512; aCDS=685; aEX=1011; %C/X=67; sC=14738558; sX=21763338 # original consensus1 gene set cacao9.consensus1 : n=34996; aCDS=1240; aEX=1869; %C/X=66; sC=43410765; sX=65414993 # est assemblies: not run thru genefindcds/pinfix ; cgbAssembly.bean : n=25502; aCDS=702; aEX=1185; %C/X=59; sC=17906790; sX=30230721 cgbAssembly.leaf : n=24590; aCDS=641; aEX=1107; %C/X=57; sC=15775000; sX=27242739 cgbAssembly.pistil1: n=20577; aCDS=470; aEX=735; %C/X=63; sC=9673284; sX=15135886 cgbAssembly.pistil2: n=26029; aCDS=593; aEX=1010; %C/X=58; sC=15460790; sX=26298700 # est asm after pinfix: cgbfix.bean : n=25502; aCDS=726; aEX=1184; %C/X=61; sC=18526730; sX=30216544 cgbfix.leaf : n=24590; aCDS=662; aEX=1107; %C/X=59; sC=16280445; sX=27224814 cgbfix.pistil1: n=20577; aCDS=496; aEX=735; %C/X=67; sC=10217475; sX=15133436 cgbfix.pistil2: n=26029; aCDS=640; aEX=1010; %C/X=63; sC=16684032; sX=26289604 #..................................... # read counts / scaffold, sc7 more of a puzzle cacao1chloroplast 160604 6204230 3863 % cacao1mito 544287 4072366 748 % cacao1rdna_scaffold_175 10639 70080890 658716 % scaffold_1 38988864 48198736 123 % scaffold_10r 25448839 17898695 70 % scaffold_2 42436413 37106286 87 % scaffold_3 34397752 35810031 104 % scaffold_4 33492547 35485292 105 % scaffold_5 40442896 36178300 89 % scaffold_6 27290986 24576787 90 % scaffold_7 24379470 14177411 58 % scaffold_8 21543242 24172042 112 % scaffold_9 42035188 46894096 111 % scaffother 15537478 13876705 89 % See read_scaf78.count for counts per read group, some groups have 1.0 ratio for 7/8, most have 0.50 .. 0.60. Chr 7 has 50% fewer rna reads than other chromosomes. Another oddity of chr7, along w/ 2x more non-align spans. Also 7 has only 70% as many introns (proportional to # reads) There may also be an EST deficit on Chr 7, which has lowest EST read count also: sc6.longfa 558854 sc7.longfa 316387 0.65 s7/s8 sc8.longfa 486875