My Blog

0. BIOPYTHON mafIO ==================     def mafParsing(self,chrName ):         idx = MafIndex(chrName + ".mafindex", chrName + ".maf", "hg19." + chrName )         multiple_alignment = idx.get_spliced([60000] ,  [60020] , strand = "-")               AlignIO.write(multiple_alignment, "hg19_All.fa", "fasta")         total_bases_tupBel1  = 0;         for seqrec in multiple_alignment:                 if seqrec.id.startswith("tupBel1"): # mm4  tupBel1  tanvir                     # don't count gaps as bases                     total_bases_tupBel1 += len( str(seqrec.seq).replace("-", "") )                            print "total of %s base align" %total_bases_tupBel1  + " %s"   %total_bases_tupBel1 1. mafsInRegion ===================================== desc :  Find multi alignment output in different species usage :  ./mafsInRegion -keepInitialGaps   -outDir test.bed ou

return the match index /  the next index if no match int findDownstreamStart( int key)     {         int imax = vectNegStart.size();         int imin = 0;         int imid = (imin + imax) / 2;                 if( key <= vectNegStart.get(0))             return 0;         if( key >= vectNegStart.get(imax-1 ))             return imax-1;                     // continue searching while [imin,imax] is not empty         while (imax >= imin)         {             /* calculate the midpoint for roughly equal partition */             imid = (imin + imax) / 2;             // determine which subarray to search             if      ( vectNegStart.get(imid) <  key)                 // change min index to search upper subarray                 imin = imid + 1;             else if (  vectNegStart.get(imid) > key )                 // change max index to search lower subarray                 imax = imid - 1;             else                 // key found at in

// create map LinkedHashMap hashmapGene = new LinkedHashMap() (); // insert element into map             if( hashmapGene.containsKey(geneID) )             {                                GeneInfo gene = hashmapGene.get(geneID);                 gene.determineClass(tmp[3]);                            }else             {                 hashmapGene.put( geneID , new GeneInfo(chrom , st,end, promStartCoord, promEndCoord ,                         tmp[3],tmp[4],strand ,tmp[6],tmp[7],tmp[8],tmp[9],tmp[10],tmp[11] ,                         strLine  ) );             }       // retrieve element from map                  Set set = hashmapGene.entrySet();             System.out.println("Total Unique entry:" + set.size() ) ;             Iterator itr = set.iterator();             while(itr.hasNext()){                 Map.Entry me = (Map.Entry) itr.next();                 String id = (String)me.getKey();                 GeneInfo geneObj = (GeneInfo) me