I see!

2013-09-10

How to install LaTeX packages

Since I tend to forget how easy it is to install LaTeX packages, here is a very short how-to:

Go to http://www.ctan.org/tex-archive/install/macros/latex/contrib and find the package you want
download the .tds.zip file
change to the folder were your downloaded file was saved to
unzip it to the directory were your texmf is stored, for me (Ubuntu) i.e.: /home/username/texmf/
update tex, in the terminal, type: texhash ~/texmf/
Done!

Or, in the terminal with the glossaries package as an example:

cd ~/Downloads/

wget http://www.ctan.org/tex-archive/install/macros/latex/contrib/glossaries.tds.zip

unzip ~/Downloads/glossaries.tds.zip -d ~/texmf/

texhash ~/texmf/

Note: This way, the new package is only installed for you, not for all users!

2013-09-05

Create a Multiple Sequence Alignment with masked parts and replace masked parts by propper sequence

The reason why I wrote the following Script is that I wanted to align a couple of sequences which contained repeats. Not surprisingly the alignment muscle produced did not look very nice, due to those repeats.

Unfortunately, one cannot tell muscle to ignore all small caps letters (afaik), hence I replaced all small caps by 'N' s. Making again an alignment produced a much better result. However, my sequences still contain the 'N's, of course. This is where the small Script kicks in: It replaces all N's with the character that is supposed to be found at this position.

Here is the code of the Script:

#! /usr/bin/python

"""
Script to "unmask" an aligned sequence.

Takes an alignment file of masked sequences and a file with the unmasked sequences
and returns an alignment file with unmasked sequences.

Written by: Tabea Hoehmann
"""

import sys
from textwrap import wrap

def read_fasta_file(fasta_filename):
"""
Read a fasta file and save the sequences to a dictionary.
Returns a dictionary with the sequence id as key and the
sequence as value.
"""
seq={}
fasta_filehandle = open(fasta_filename, 'r')
i=0
for line in fasta_filehandle.readlines():
if '>' in line:
name = line.strip().replace('>','').upper() ### !!!! Only the first element !!!!
seq[name] = ""
else:
seq[name] += line.strip()
if not set(line.strip()).issubset(set('actgkrynACTGKRYN-')):
print("Unknown character(s) in line %s of sequence file %s: %s\nScript proceeds as if nothing happened..."%(i, fasta_filename, ", ".join(set(line)-set('actgkrynACTGKRYN-'))))
i+=1
fasta_filehandle.close()
return(seq)

# Read fasta files
try:
masked_aln_fasta = read_fasta_file(sys.argv[1])
unmasked_fasta = read_fasta_file(sys.argv[2])
except IndexError:
print("Usage: unmask_fasta_aln.py file1 file 2")
print("file1: The alignment file that is masked by N's")
print("file1: The fasta file that is not masked")
sys.exit()

# New dictionary to store edited sequences
new_unmasked_aln_fasta = {}

# Iterate over every sequence, replace all N's by nucleotides
for masked_id, masked_seq in masked_aln_fasta.items():
unmasked_seq = unmasked_fasta[masked_id]
new_unmasked_aln_seq = ""
nt_pos = 0
for nt in masked_seq:
if nt not in 'nN':
new_unmasked_aln_seq += nt
else:
new_unmasked_aln_seq += unmasked_seq[nt_pos]
if nt != '-':
nt_pos+=1
new_unmasked_aln_fasta[masked_id] = new_unmasked_aln_seq

# Print result to stdout
for id, seq in new_unmasked_aln_fasta.items():
print('>%s\n%s'%(id, '\n'.join(wrap(seq, 100, break_on_hyphens=False))))

2013-07-29

Use cut to extract a certain text range from a one-line DNA sequence

I wanted to get the characters at position 27947 to 28203 from a text file containing just one DNA sequence as a very long line. This is the cut command that worked:

cut -c 27947-28203 file.txt

2013-07-19

How to automatically create Venn diagram for five sets

Wikipedia has a very nice Venn diagram for five sets. I use this svg file to automatically create a Venn diagram displaying the number of intersecting elements of five sets.

Number of intersecting elements

For this we need the number of intersecting elements for all possible set combinations, in my example it looks like this:

numbers.txt

A 13644

B 14729

C 14690

D 13725

E 13742

AB 3689

AC 3616

AD 3523

AE 3496

BC 14281

BD 12852

BE 12694

CD 13215

CE 13060

DE 13563

ABC 3609

ABD 3507

ABE 3480

ACD 3513

ACE 3487

ADE 3496

BCD 12849

BCE 12694

BDE 12694

CDE 13056

ABCD 3506

ABCE 3480

ABDE 3480

ACDE 3487

BCDE 12694

ABCDE 3480

To get these numbers I used a small R script and formatted the output with a little bit of TextWrangler and grep to look like in the file numbers.txt

Here is the short R script:

A <- read.delim('Info1.csv')$name

B <- read.delim('Info2.csv')$name

C <- read.delim('Info3.csv')$name

D <- read.delim('Info4.csv')$name

E <- read.delim('Info5.csv')$name

# all sets

print(paste("A",length(A)))

print(paste("B",length(B)))

print(paste("C",length(C)))

print(paste("D",length(D)))

print(paste("E",length(E)))

# all combinations of two sets

print(paste("A, B", length(intersect(A, B))))

print(paste("A, C", length(intersect(A, C))))

print(paste("A, D", length(intersect(A, D))))

print(paste("A, E", length(intersect(A, E))))

print(paste("B, C", length(intersect(B, C))))

print(paste("B, D", length(intersect(B, D))))

print(paste("B, E", length(intersect(B, E))))

print(paste("C, D", length(intersect(C, D))))

print(paste("C, E", length(intersect(C, E))))

print(paste("D, E", length(intersect(D, E))))

# all combinations of three sets

print(paste("A, B, C", length(intersect(A, intersect(B, C)))))

print(paste("A, B, D", length(intersect(A, intersect(B, D)))))

print(paste("A, B, E", length(intersect(A, intersect(B, E)))))

print(paste("A, C, D", length(intersect(A, intersect(C, D)))))

print(paste("A, C, E", length(intersect(A, intersect(C, E)))))

print(paste("A, D, E", length(intersect(A, intersect(D, E)))))

print(paste("B, C, D", length(intersect(B, intersect(C, D)))))

print(paste("B, C, E", length(intersect(B, intersect(C, E)))))

print(paste("B, D, E", length(intersect(B, intersect(D, E)))))

print(paste("C, D, E", length(intersect(C, intersect(D, E)))))

# all combinations of four sets

print(paste("A, B, C, D", length(intersect(A, intersect(B, intersect(C, D))))))

print(paste("A, B, C, E", length(intersect(A, intersect(B, intersect(C, E))))))

print(paste("A, B, D, E", length(intersect(A, intersect(B, intersect(D, E))))))

print(paste("A, C, D, E", length(intersect(A, intersect(C, intersect(D, E))))))

print(paste("B, C, D, E", length(intersect(B, intersect(C, intersect(D, E))))))

# combination of five sets

print(paste("ABCDE", length(intersect(A, intersect(B, intersect(C, intersect(D, E)))))))

Getting these numbers into the diagram

Since the svg file only consists of 51 lines, I display it here:

<?xml version="1.0" encoding="UTF-8"?>

<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">

<title>Radially-symmetrical Five-set Venn Diagram</title>

<desc>Devised by Branko Gruenbaum and rendered by CMG Lee.</desc>

<defs>

</g>

</defs>

<text x="0" y="0" dy="0.7ex">ABCDE</text>

</g>

</svg>

You can see, that the names of the different overlapping fields are annotated as "A", "B", "AB", and so forth. Hence, I can just replace them with sed, to display my own data!

This I do with the following bash/sed command:

while read LINE

NAME=$( echo $LINE | awk '{print $1 }')

NUMB=$( echo $LINE | awk '{print $2 }')

echo $NAME $NUMB

sed -i tmp "s/\>${NAME}\</\>${NUMB}\</g" venn.svg

echo "sed -i tmp "s/\>\${NAME}\</\>\${NUMB}\</g" venn.svg"

done < 'numbers.txt'

The altered Venn diagram

The output is indeed exactly what I wanted:

Of course, I still need to adapt the font size, but thanks to sed, this won't be a big issue.

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 see stdout output in the terminal and pipe it into a file

Since I always have to look up how to print the stdout output to the screen AND to a file at the same time, I put this small command here:

command 2>&1 | tee outputfile.txt

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))

2013-02-18

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

2013-01-30

RepeatMasker ambiguous bases

When trying to use a chromosome fasta file with RepeatMasker, one might gets the error message that "ambiguous bases" were detected because the fasta file contains long stretches of N's.
In order to avoid this problem, you can replace the "N" characters by "X"'s before starting RepeatMasker.

2012-11-13

How to split a string in R using strsplit

Since the first google and Rseek.org results did not help me with my problem of splitting a string in R, I want to clip it here, so I do not have to click through all of these pages again :)

The command is as easy as this:

> y <- "100,200,300"
> x <- strsplit(y, ',')
> x
[[1]]
[1] "100" "200" "300"

2012-10-31

R: Function to perform t-test on all possible combinations of columns in a dataframe

Sometimes I need to perform a t-test on all possible combinations of the columns of a data-frame and rather than computing them all one by one, I wrote this small function:

ttest_all_possible <- function(dataframe) {
    nr_columns = dim(dataframe)[2]
    results <- matrix(NA, nr_columns, nr_columns)
    for (i in seq(1, nr_columns)){
        for (j in seq(1, nr_columns)){
            results[i, j] <- t.test(dataframe[,i], dataframe[,j])$p.value
        }
    }
    return(results)
}

Works fine:
> test <- cbind(rnorm(10),rnorm(10), rnorm(10))
> test
             [,1]       [,2]       [,3]
[1,] 0.64061950 -0.5643179 0.7450485
[2,] 1.90703328 -0.8000197 0.2444611
[3,] 0.96346636 -0.6278075 0.1283072
[4,] -0.25758024 -0.9657888 -0.8832595
[5,] 0.36934381 1.4149775 0.8837320
[6,] 0.60584130 0.2556500 -0.2437533
[7,] -0.24559562 -0.7774171 0.1279253
[8,] -1.53276425 0.3557177 0.1611155
[9,] -1.55474156 0.3488891 0.1025763
[10,] -0.08366046 -0.2823561 1.3618903
> ttest_all_possible(test)
          [,1]      [,2]      [,3]
[1,] 1.0000000 0.5603004 0.6495045
[2,] 0.5603004 1.0000000 0.1817692
[3,] 0.6495045 0.1817692 1.0000000

2012-10-30

How to substutitute characters in all files of a directory using regex and perl

perl -i -pe 's/CharacterToSubstiture/SubstituteWithThis/g' myfilename

Example:
Imagine I have got these two files:

shopping_list_groceries.txt:
- 2 apples
- 3 red paprika

shopping_list_clothes.txt:

- 2 pairs of socks

- 1 red t-shirt

And now I decide to not like red any more, but to prefer green:

perl -i -pe 's/red/green/g' shopping_list*txt

The text files now look like this:
shopping_list_groceries.txt:
- 2 apples
- 3 green paprika

shopping_list_clothes.txt:

- 2 pairs of socks

- 1 green t-shirt

What does the one-liner do?

perl # use perl :)

-i # change in the file, if you skip this part, the changes will be printed to stdout

-pe # I don't know

's/red/green/g' # RegEx:

's # Substitute...

/red # ...every occurence of the string "red"...

/green # ...with the string "green"...

/g' # ... everywhere ...

shopping_list*.txt # ...in all files starting with 'shopping_list' and ending with '.txt

2012-10-02

How to print all the attributes of an object in python

from pprint import pprint
pprint (vars(my_object))     # e.g. class instance

This prints all the attributes/properties this object currently has in a neat way.

2012-09-19

How to create a tar archive

tar cvzf archive.tar.gz *files*

2012-09-18

Bash Shell: Check if a file exists, and has a size larger than 0

To check, whether a file exists and has a file size larger than 0 use this expression:

if [ -s someFile.txt ]
then
...
fi

Bash Shell: Check if a file exists

To check, whether a file exists use this expression:

First, create the file thisFileExists.txt:
touch thisFileExists.txt

Now, check if it exists:
if [ -f thisFileExists.txt ]
then
echo It is there
fi

This is what it looks like in terminal:

Bash Shell: Repeat a command every x seconds, e.g. to update list of files

When running a command that changes/creates files and takes a while, one might be interested in seeing how progress is going.
For this the "watch" command can be useful:

watch -n 2 ls *file

This will print the output of "ls *file" every 2 seconds.

watch -n 2 -d ls *file

Adding -d will highlight all the changes in comparison to the last time the command was performed

2012-07-28

Creating an animated gif image

Imagemagick can also be used to create animations!

convert -delay 10 -loop 0 DSC003* animation.gif

Where DSC003* are the files that are going to be used for the animations. They can be in about any format.

As a (rather stupid) example, I created 16 barplots in R with different colors and animated this:

in R:
$ pdf('example.pdf')
$ for (i in seq(5,20)){barplot(1:20,col=rainbow(i))}
$ dev.off()

in Shell:
$ convert example.pdf example.png
$ convert -delay 10 -loop 1 example-*png animation.gif

The 16 different barplots.

The slightly psychedelic animation.

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

2012-07-17

How to get nice Python Syntax Highlighting in LaTeX documents

Pygmentize is a command that can be performed on a python Script to produce a .tex document with nice syntax highlighting

pygmentize -f latex script.py > script.tex

2012-07-14

How to install .deb packages on Linux

$ sudo dpkg -i package.deb

2012-07-08

Pickle module in Python

The pickle module in Python 2.3 can be used to save a dictionary in a file and makes it accessible for later.

import pickle

((I will add more information later))