How to convert a BED file with one line per gene/item into a BED file with one line per block

If a BED file  looks for example like this

chr1    11873    14361    AM992877    0    +    11873    14361    0    3    354,109,1141,    0,739,1347,
chr1    11873    14361    AM992881    0    +    11873    14361    0    3    354,127,959,    0,721,1529,
chr1    11873    14362    AM992878    0    +    11873    14362    0    5    354,52,436,299,404,    0,772,1347,1785,2085,
chr1    11873    14362    AM992879    0    +    11873    14362    0    3    354,109,1142,    0,739,1347,
chr1    11873    14409    AM992871    0    +    11873    14409    0    3    354,109,1189,    0,739,1347,

it contains one line for every mRNA (in this case), giving the mRNA start and end coordinates in the first three columns. However, the blocks' coordinates are given in the last two columns and if you use a tool such as, e.g., Bedtools, you can only find results for the whole length of the mRNA.
I wanted to get results based on the individual blocks of the mRNA.

To get these I wrote a script that turns the above BED file into this format:
chr1    11873    12227    AM992877    0    +    11873    14361    0    3    354,109,1141,    0,739,1347,
chr1    12612    12721    AM992877    0    +    11873    14361    0    3    354,109,1141,    0,739,1347,
chr1    13220    14361    AM992877    0    +    11873    14361    0    3    354,109,1141,    0,739,1347,
chr1    11873    12227    AM992881    0    +    11873    14361    0    3    354,127,959,    0,721,1529,
chr1    12594    12721    AM992881    0    +    11873    14361    0    3    354,127,959,    0,721,1529,
chr1    13402    14361    AM992881    0    +    11873    14361    0    3    354,127,959,    0,721,1529,
chr1    11873    12227    AM992878    0    +    11873    14362    0    5    354,52,436,299,404,    0,772,1347,1785,2085,
chr1    12645    12697    AM992878    0    +    11873    14362    0    5    354,52,436,299,404,    0,772,1347,1785,2085,
chr1    13220    13656    AM992878    0    +    11873    14362    0    5    354,52,436,299,404,    0,772,1347,1785,2085,
chr1    13658    13957    AM992878    0    +    11873    14362    0    5    354,52,436,299,404,    0,772,1347,1785,2085,
chr1    13958    14362    AM992878    0    +    11873    14362    0    5    354,52,436,299,404,    0,772,1347,1785,2085,
chr1    11873    12227    AM992879    0    +    11873    14362    0    3    354,109,1142,    0,739,1347,
chr1    12612    12721    AM992879    0    +    11873    14362    0    3    354,109,1142,    0,739,1347,
chr1    13220    14362    AM992879    0    +    11873    14362    0    3    354,109,1142,    0,739,1347,
chr1    11873    12227    AM992871    0    +    11873    14409    0    3    354,109,1189,    0,739,1347,
chr1    12612    12721    AM992871    0    +    11873    14409    0    3    354,109,1189,    0,739,1347,
chr1    13220    14409    AM992871    0    +    11873    14409    0    3    354,109,1189,    0,739,1347,

This is the script, which works nice, if the BED file has the exact same format as in the example:


import sys

infile = sys.argv[1]

infile_handle = open(infile)

for line in infile_handle.readlines():
    chr, chr_start, chr_end, name, score, strand, thickStart, \
    thickEnd, itemRGB, blockCount, blockSizes, blockStarts = \
    line.split('\t')
    chr_start = int(chr_start)
    chr_end   = int(chr_end)
    for i, blockSize in enumerate(blockSizes.split(',')[:-1]):
        blockSize = int(blockSize)
        blockStart = int(blockStarts.split(',')[:-1][i])
        s = "\t".join([chr, str(chr_start + blockStart), str(chr_start + blockStart + blockSize), name, score, strand, thickStart, thickEnd, itemRGB, blockCount, blockSizes, blockStarts])
        print s,

How to convert fasta file to uppercase/lowercase, without altering the annotation line

This awk-oneliner converts all lines which do not contain a '>' in the beginning to upper case:

$ awk '{ if ($0 !~ />/) {print toupper($0)} else {print $0} }' name.fasta

How to convert .maf alignment format to a RepeatMasker-output-file-like format

The following python script takes a .maf alignment file, which looks for example like this:


# lastz.v1.02.00 --format=maf --output=NM_144963.2-3.maf 
#
# hsp_threshold      = 3000
# gapped_threshold   = 3000
# x_drop             = 910
# y_drop             = 9400
# gap_open_penalty   = 400
# gap_extend_penalty = 30
#        A    C    G    T
#   A   91 -114  -31 -123
#   C -114  100 -125  -31
#   G  -31 -125  100 -114
#   T -123  -31 -114   91
a score=11361
s target 0 133 +  5344 GGCGGGAGGCTTAGGGTGGGTCGCGGTCGCGGCCCCCGCCGCAGGCCCTGGCGGCGGCATCGTGAACATAGACGTGGAGTTCCACATCCGGCACAACTACCCCTGGAACAAGTTGCCGGCCAACGTGAGGCAG
s query           0 133 + 46809 GGCGGGAGGCGTAGTGTGGGTCGCGGTGGCGGCCCCGGCCGCCGGCCCTGGCCGCGGCATCATGAACATAGACGTGGAGTTCCACATCCGGCACAACTACCCCTGGAACAAGTTGCCGGCCAACGTGAGACAG

a score=7552
s target  207 151 +  5344 tacaggaagtattggtcttgctgttcagagtggcagaggagaaaga--------aaatcccagcC---GTCTTGGAAATTCACAGAGAGAATATGAAAATCAGGTTGTCCTGTACAGTATCTGCAATCAGTTACGATACAGAAATAACTTAAACATGTCAAG
s query           5371 162 + 46809 tatgggaaaaaatagtactattgttcagtttggaagttgtttaaaatctttttgaaatttctcctagAGTCTTGGAAATTCACAGAGAGAATATGAAAAGCAGGTTGTCCTGTACAGTATCCGCAATCAGTTACGATATAGAAATAACTTAGgtaagtagag

a score=13285
s target  344 152 +  5344 TTAAACATGTCAAGAAAGATGAACGCAGATACTATGAGGAACTGCTAAACTACAGCCGAGATCATCTCATGCTGTACCCTTACCGTCTATTGGATATTATGGTGAAGAGCTTGAGCATAACACCATTTTCATATTACACTGGGATTATGGAG
s query           6505 152 + 46809 TTAAACATGTCAAGAAAGATGAACGCAGATACTATGAGGAACTGCTAAAGTACAGCCGAGATCATCTCATGCTGTACCCTTACCATCTATCGGATATTATGGTGAAAGGCTTGAGGATAACACCATTTTCATATTATACTGGGATTATGGAG

a score=6190
s target  536 68 +  5344 GTCTAAGGCTTCTTGGCATAGGAAGAAACCAGTATATTGATCTTATGATTCAGTGTAGATCATCAAAA
s query           9383 68 + 46809 GTCTAAGGCTTCTTGGCATAGGAAGAAACCAGTATATTGATCTTATGAATCAGTGTAGATCATCAAAA


and converts it into a RepeatMasker-output-like format that looks like this:
SW      subst   del_in_retro    insert_in_retro query   qstart  qend    qleft   compl   match   match_class     mstart  mend    mleft   id
11361 nA nA nA hg19_ucscRetroAli5_NM_144963.2-3 0 133 5211 +,+ hg19_refGene_NM_144963 nA 0 133 46676 1
7552 nA nA nA hg19_ucscRetroAli5_NM_144963.2-3 207 358 4986 +,+ hg19_refGene_NM_144963 nA 5371 5533 41276 2
13285 nA nA nA hg19_ucscRetroAli5_NM_144963.2-3 344 496 4848 +,+ hg19_refGene_NM_144963 nA 6505 6657 40152 3
6190 nA nA nA hg19_ucscRetroAli5_NM_144963.2-3 536 604 4740 +,+ hg19_refGene_NM_144963 nA 9383 9451 37358 4

I save also the output that will not go into the new file into variables, in case I maybe want to include it as well.




import sys
INFILENAME=sys.argv[1]

WITH_HEADER=False
try:
WITH_HEADER=sys.argv[2]
except:
pass

INFILE = open(INFILENAME, 'r')

ALNS=[]
i=0
for LINE in INFILE.readlines():
if 'a score=' in LINE:
SCORE=LINE.split('=')[1].strip()
i+=1
QUERY_PASSED=False
MATCH_PASSED=False
elif LINE.startswith('s') and not QUERY_PASSED:
   MUERY, MSTART, MSIZE, MSTRAND, MSRCSIZE, MSEQUENCE = LINE.split()[1:]
   MEND = str(int(MSTART) + int(MSIZE))
   MLEFT = str(int(MSRCSIZE) - int(MEND))
   MUERY_PASSED = True
elif LINE.startswith('s') and QUERY_PASSED:
   QATCH, QSTART, QSIZE, QSTRAND, QSRCSIZE, QSEQUENCE = LINE.split()[1:]
   QEND = str(int(QSTART) + int(QSIZE))
   QLEFT = str(int(QSRCSIZE) - int(QEND))
   QATCH_PASSED=True
   ALNS.append((SCORE, 'nA', 'nA', 'nA', QUERY, QSTART, QEND, QLEFT, MSTRAND, MATCH, 'nA', MSTART, MEND, MLEFT, str(i)))

INFILE.close()

if WITH_HEADER:
print("SW      subst   del_in_retro    insert_in_retro query   qstart  qend    qleft   compl   match   match_class     mstart  mend    mleft   id")    
for LINE in ALNS:
print('\t'.join(LINE))
 

Use Imagemagick to swap blue and red colors

For a presentation I wanted an image to appear in magenta, rather than in cyan. This is how the image originally looks like:




blue.jgp

With the following command I separated the R, G, B color channels and swapped the Red and the Green channel, to produce a picture in which all blue colors appear in red:

$ convert blue.jpg -separate -swap 0,1 -combine red.jpg

 And that is how the result looks like:




red.jpg


I got this information from: http://www.imagemagick.org/Usage/color_basics/#combine

How to bulk resize JPEG pictures via Imagemagick

Sometimes you need to reduce the size of many pictures, for example to share them on the internet.

A nice way of doing this is via Shell:

convert FILE.JPG -resize 2448x1836 ./smaller_pictures/FILE.JPG

Or by using a for loop:

for F in *.JPG

do

convert $F -resize 2448x1836 ./smaller_pictures/$F
done