2013年4月2日星期二

modify sequence name in a fasta file

cat genes.fasta | awk '{if (/^>/) {print $1} else {print $0}}' > genes2.fasta

To shorten the sequence names in fasta file.

change
>sequense1 human person1
ATGC

to

>sequense1
ATGC

Some unix/perl oneliners for Bioinformatics

http://genomics-array.blogspot.ca/2010/11/some-unixperl-oneliners-for.html

   File format conversion/line counting/counting number of files etc.

1.    $ wc –l   : count number of lines in a file.
2.    $ ls | wc –l        : count number of files in a directory.
3.    $ tac     : print the file in reverse order e.g; last line first, first line last.
4.    $ rev     : reverse the file in lines.
5.    $ sed 's/.$//' or sed 's/^M$//' or sed 's/\x0D$//' : converts a dos file into unix mode.
6.    $sed "s/$/`echo -e \\\r`/" or sed 's/$/\r/' or sed "s/$//": converts a unix newline into a DOS newline.
7.    $ awk '1; { print "" }' : Double space a file.
8.    $ awk '{ total = total + NF }; END { print total+0 }' : prints the number of words in a file.
9.    $sed '/^$/d' or [grep ‘.’] : Delete all blank lines in a file.
10.    $sed '/./,$!d' : Delete all blank lines in the beginning of the file.
11.    $sed -e :a -e '/^\n*$/{$d;N;ba' -e '}': Delete all blank lines at the end of the file.
12.    $sed -e :a -e 's/<[^>]*>//g;/
13.    $sed 's/^[ \t]*//' : deleting all leading white space tabs in a file.
14.    $ sed 's/[ \t]*$//' : Delete all trailing white space and tab in a file.
15.    $ sed 's/^[ \t]*//;s/[ \t]*$//' : Delete both leading and trailing white space and tab in a file.


2.2    Working with Patterns/numbers in a sequence file
16.    $awk '/Pattern/ { n++ }; END { print n+0 }' : print the total number of lines containing the word pattern.
17.    $sed 10q : print first 10 lines.
18.    $sed -n '/regexp/p' : Print the line that matches the pattern.
19.    $sed '/regexp/d' : Deletes the lines that matches the regexp.
20.    $sed -n '/regexp/!p' : Print the lines that does not match the pattern.
21.    $sed '/regexp/!d' : Deletes the lines that does NOT match the regular expression.
22.    $sed -n '/^.\{65\}/p' : print lines that are longer than 65 characters.
23.    $sed -n '/^.\{65\}/!p' : print lines that are lesser than 65 characters.
24.    $sed -n '/regexp/{g;1!p;};h' : print one line before the pattern match.
25.    $sed -n '/regexp/{n;p;}' : print one line after the pattern match.
26.    $sed -n '/^.\{65\}/ {g;1!p;};h' < sojae_seq > tmp : print the names of the sequences that are larger than 65 nucleotide long.
27.    $sed -n '/regexp/,$p' : Print regular expression to the end of file.
28.    $sed -n '8,12p' : print line 8 to 12(inclusive)
29.    $sed -n '52p' : print only line number 52.
30.    $seq ‘/pattern1/,/pattern2/d’ < inputfile > outfile : will delete all the lines between pattern1 and pattern2.
31.    $sed ‘/20,30/d’ < inputfile > outfile : will delete all lines between 20 and 30.   OR sed ‘/20,30/d’ < input > output will delete lines between 20 and 30.
32.    awk '/baz/ { gsub(/foo/, "bar") }; { print }' : Substitute foo with bar in lines that contains ‘baz’.
33.    awk '!/baz/ { gsub(/foo/, "bar") }; { print }' : Substitute foo with bar in lines that does not contain ‘baz’.
34.    grep –i –B 1 ‘pattern’ filename > out : Will print the name of the sequence and the sequence having the pattern in a case insensitive way(make sure the sequence name and the sequence each occupy a single line).
35.    grep –i –A 1 ‘seqname’ filename > out : will print the sequence name as well as the sequence into file ‘out’.  

2.3    Inserting Data into a file:

36.    gawk --re-interval 'BEGIN{ while(a++<49 1="" amp="" filename="" nbsp="" s="" sub="" x="">> fileout : will insert 49 ‘X’ in the sixth position of every line.


37.    gawk --re-interval 'BEGIN{ s="YourName" }; { sub(/^.{6}/,"&" s) }; 1' : Insert your name at the 6 th position in every line.

3.    Working with Data Files[Tab delimited files]:

3.1    Error Checking and data handling:
38.    awk '{ print NF ":" $0 } ' : print the number of fields of each line followed by the line.
39.    awk '{ print $NF }' : print the last field of each line.
40.    awk 'NF > n' : print every line with more than ‘n’ fields.
41.    awk '$NF > n' : print every line where the last field is greater than n.
42.    awk '{ print $2, $1 }' : prints just first 2 fields of a data file in reverse order.
43.    awk '{ temp = $1; $1 = $2; $2 = temp; print }' : prints all the fields in the correct order except the first 2 fields.
44.    awk '{ for (i=NF; i>0; i--) printf("%s ", $i); printf ("\n") }' : prints all the fields in reverse order.
45.    awk '{ $2 = ""; print }' : deletes the 2nd field in each line.
46.    awk '$5 == "abc123"' : print each line where the 5th field is equal to ‘abc123’.
47.    awk '$5 != "abc123"' : print each line where 5th field is NOT equal to abc123.
48.    awk '$7  ~ /^[a-f]/' : Print each line whose 7th field matches the regular expression.
49.    awk '$7 !~ /^[a-f]/' : print each line whose 7th field does NOT match the regular expression.
50.    cut –f n1,n2,n3.. > output file : will cut n1,n2,n3 columns(fields) from input file and print the output in output file. If delimiter is other than TAB then give additional argument such as cut –d ‘,’ –f n1,n2.. inputfile > out
51.    sort –n –k 2,2 –k 4,4 file > fileout : Will conduct a numerical sort of column 2, and then column 4. If –n is not specified, then, sort will do a lexicographical sort(of the ascii value).

4.    Miscellaneous:
52.    uniq –u inputfile > out : will print only the uniq lines present in the sorted input file.
53.    uniq –d inputfile > out : will print only the lines that are in doubles from the sorted input file.
54.    cat file1 file2 file3 … fileN > outfile : Will concatenate files back to back in outfile.
55.    paste file1 file2 > outfile : will merge two files horizontally. This function is good for merging with same number of rows but different column width.
56.    !:p : will print the previous command run with the ‘pattern’ in it.
57.    !! : repeat the last command entered at the shell.
58.    ~ : Go back to home directory
59.     echo {a,t,g,c}{a,t,g,c}{a,t,g,c}{a,t,g,c} : will generate all tetramers using ‘atgc’. If you want pentamers/hexamers etc. then just increase the number of bracketed entities.NOTE: This is not a efficient sequence shuffler. If you wish to generate longer sequences then use other means.
60.    kill -HUP ` ps -aef | grep -i firefox | sort -k 2 -r | sed 1d | awk ' { print $2 } ' ` : Kills a hanging firefox process.
61.    csplit -n 7 input.fasta '/>/' '{*}' : will split the file ‘input.fasta’ wherever it encounters delimiter ‘>’. The file names will appear as 7 digit long strings.
62.    find . -name data.txt –print: finds and prints the path for file data.txt.
Sample Script to make set operations on sequence files:
63.    grep ‘>’ filenameA > list1  # Will list just the sequence names in a file names.
grep ‘>’ filenameB > list2 # Will list names for file 2
cat list1 list2 > tmp # concatenates list1 and list2 into tmp
sort tmp > tmp1 # File sorted
uniq –u tmp1 > uniq    # AUB – A ∩ B (OR (A-B) U (B-A)) 
uniq –d tmp1 > double  # Is the intersection (A ∩ B)
cat uniq double > Union # AUB
cat list1 double > tmp
sort tmp | uniq –u > list1uniq # A - B
cat list2 double > tmp
sort tmp | uniq –u > list2uniq # B - A     



PERL ONELINERS:

1.    perl -pe '$\="\n"'   : double space a file
2.    perl -pe '$_ .= "\n" unless /^$/' : double space a file except blank lines
3.    perl -pe '$_.="\n"x7' : 7 space in a line.
4.    perl -ne 'print unless /^$/' : remove all blank lines
5.    perl -lne 'print if length($_) < 20' : print all lines with length less than 20.
6.    perl -00 -pe '' : If there are multiple spaces, delete all leaving one(make the file a single spaced file).
7.    perl -00 -pe '$_.="\n"x4' : Expand single blank lines into 4 consecutive blank lines
8.    perl -pe '$_ = "$. $_"': Number all lines in a file
9.    perl -pe '$_ = ++$a." $_" if /./' : Number only non-empty lines in a file
10.    perl -ne 'print ++$a." $_" if /./' : Number and print only non-empty lines in a file
11.    perl -pe '$_ = ++$a." $_" if /regex/' ; Number only lines that match a pattern
12.    perl -ne 'print ++$a." $_" if /regex/' : Number and print only lines that match a pattern
13.    perl -ne 'printf "%-5d %s", $., $_ if /regex/' : Left align lines with 5 white spaces if matches a pattern (perl -ne 'printf "%-5d %s", $., $_' : for all the lines)
14.    perl -le 'print scalar(grep{/./}<>)' : prints the total number of non-empty lines in a file
15.    perl -lne '$a++ if /regex/; END {print $a+0}' : print the total number of lines that matches the pattern
16.    perl -alne 'print scalar @F' : print the total number fields(words) in each line.
17.    perl -alne '$t += @F; END { print $t}' : Find total number of words in the file
18.    perl -alne 'map { /regex/ && $t++ } @F; END { print $t }' : find total number of fields that match the pattern
19.    perl -lne '/regex/ && $t++; END { print $t }' : Find total number of lines that match a pattern
20.    perl -le '$n = 20; $m = 35; ($m,$n) = ($n,$m%$n) while $n; print $m' : will calculate the GCD of two numbers.
21.    perl -le '$a = $n = 20; $b = $m = 35; ($m,$n) = ($n,$m%$n) while $n; print $a*$b/$m' : will calculate lcd of 20 and 35.
22.    perl -le '$n=10; $min=5; $max=15; $, = " "; print map { int(rand($max-$min))+$min } 1..$n' : Generates 10 random numbers between 5 and 15.
23.    perl -le 'print map { ("a".."z",”0”..”9”)[rand 36] } 1..8': Generates a 8 character password from a to z and number 0 – 9.
24.    perl -le 'print map { ("a",”t”,”g”,”c”)[rand 4] } 1..20': Generates a 20 nucleotide long random residue.
25.    perl -le 'print "a"x50': generate a string of ‘x’ 50 character long
26.    perl -le 'print join ", ", map { ord } split //, "hello world"': Will print the ascii value of the string hello world.
27.    perl -le '@ascii = (99, 111, 100, 105, 110, 103); print pack("C*", @ascii)': converts ascii values into character strings.
28.    perl -le '@odd = grep {$_ % 2 == 1} 1..100; print "@odd"': Generates an array of odd numbers.
29.    perl -le '@even = grep {$_ % 2 == 0} 1..100; print "@even"': Generate an array of even numbers
30.    perl -lpe 'y/A-Za-z/N-ZA-Mn-za-m/' file: Convert the entire file into 13 characters offset(ROT13)
31.    perl -nle 'print uc' : Convert all text to uppercase:
32.    perl -nle 'print lc' : Convert text to lowercase:
33.    perl -nle 'print ucfirst lc' : Convert only first letter of first word to uppercas
34.    perl -ple 'y/A-Za-z/a-zA-Z/' : Convert upper case to lower case and vice versa
35.    perl -ple 's/(\w+)/\u$1/g' : Camel Casing
36.    perl -pe 's|\n|\r\n|' : Convert unix new lines into DOS new lines:
37.    perl -pe 's|\r\n|\n|' : Convert DOS newlines into unix new line
38.    perl -pe 's|\n|\r|' : Convert unix newlines into MAC newlines:
39.    perl -pe '/regexp/ && s/foo/bar/' : Substitute a foo with a bar in a line with a regexp.

Some other Perl Tricks

Want to display some progress bars while perl does your job:

For this perl provides a nice utility called "pipe opens" ('perldoc -f open' will provide more info)
open(my $file, '-|', 'command','option', 'option', ...) or die "Could not run tar ... - $!";
  while (<$file>) {
       print "-";
  }
  print "\n";
  close($file);
 
Will print - on the screen till the process is completed 

creating password less connections to multiple remote machines


creating password less connections to multiple remote machines


Many times while working in multiple machines, you may like to automate certain processes where the programs can directly access information/data from another server effortlessly - without a password. This can be done using ssh-keygenprotocol.

So, what happens here is; you have a local machine, lets call it 'A' and you have a remote machine, lets call it as 'B'. You have an account in 'B' and that is say 'myname'. Everytime you log into that machine using ssh, you have to do something like:
$ ssh myname@B
$Password:
$myname@B:

In order to directly log into a machine without a password, you have to generate a pair of keys; called as a public key and a private key. The public key is the public information and the private key is only known to your local machine i.e; 'A'. You can use ssh-keygen to create a pair of keys in a given time. This is how you should proceed:
$ ssh-keygen -t rsa
Generating public/private rsa key pair.
Enter file in which to save the key (/home/Sucheta/.ssh/id_rsa): /home/Sucheta/.
ssh/iicb_rsa [ Remember here to enter a new file name, else it will over write any other id_rsa file that you may have saved earlier for any other computer ]
Enter passphrase (empty for no passphrase): [Enter a paraphrase that is > 4 characters long. This is essential if your computer has more than one pair of public private keys for more than one remote server]
Enter same passphrase again:
Your identification has been saved in /home/Sucheta/.ssh/iicb_rsa.
Your public key has been saved in /home/Sucheta/.ssh/iicb_rsa.pub.
The key fingerprint is:
95:13:96:1b:66:ef:36:74:25:76:05:23:64:58:bb:94 Sucheta@Sucheta-PC
The key's randomart image is:
+--[ RSA 2048]----+
|          o== o.o|
|         .*+ +o.o|
|         o++E. + |
|         ..oo..  |
|        S  o..   |
|            +    |
|           . .   |
|                 |
|                 |
+-----------------+
Then do:
$ssh-copy-id myname@B
[This command will append your public key in the ~/.ssh/authorized_keys file in the remote host. You can also do this manually by logging back to your remote computer and copy pasting your public key in the 'authorized_keys' file. Make sure your public key is copy pasted in one single line.
Another thing to remember is, depending on the OS and version, the file that needs to have the public key in the remote machine may be different. In order to confirm that it is indeed called as "authorized_keys" do the following: 
[root@Apala ssh]# cat /etc/ssh/sshd_config | grep Keys
# HostKeys for protocol version 2
#AuthorizedKeysFile     .ssh/authorized_keys
#AuthorizedKeysCommand none
#AuthorizedKeysCommandRunAs nobody

This tells you indeed the file that stores public key in your remote computer is named as authorized_keys .

Next time you try to create another passwordless connection to another computer, just repeat the above steps. And always remember to write the public and private key into different files, else it will overwrite contents into id_rsa and id_rsa.pub file. Add a passphrase too.

One more important thing to remember is to check the file permission for "authorized_keys" file in the remote machine.Always set it to 700.

Using this, you can also automate file transfer by an sftp or any other remote ftp protocols

2013年3月31日星期日

PGML - plant genome mapping lab

http://www.plantgenome.uga.edu/personnel.html


MULTIPLE COLLINEARITY SCAN - MCSCAN


MCScanX-transposed is a software package able to detect transposed gene duplications that occurred within different epochs based on applying MCScanX within and between related genomes, also useful for integrative analysis of gene duplication modes and annotating a gene family of interest with gene duplication modes.

MCScan is an algorithm to scan multiple genomes or subgenomes to identify putative homologous chromosomal regions, then align these regions using genes as anchors. MCScanXtoolkit implements an adjusted MCScan algorithm for detection of synteny and collinearity and extends the software by incorporating 15 utility programs for display and further analyses. Compared with MCScan version 0.8, MCScanX has the following new features:

2013年3月29日星期五

The Santos Lab


>> mapping_pe_reads_w_bwa_bowtie2.sh - Shell script for mapping Illumina reads to scaffold(s) in FASTA format. Needs working installation of BWASAMtools and Bowtie 2.
>> MPI-enabled MrBayes and PhyML - Precompiled binaries of the phylogenetic programs MrBayes and PhyML capable of utilizing multiple CPUs simultaneously (built for Apple Intel systems).
>> remote_blast_client.prl - Performs various BLAST searches against NCBI's databases.
>> blast_parse_all.prl - Parses BLAST reports for all HSPs with BioPerl's Bio::SearchIO module.
>> blast_parse_single.prl - Parses BLAST reports for single best HSP with BioPerl's Bio::SearchIO module.
>> blast2ps.prl - Creates a graphical representation of BLAST reports as a Postscript file.
>> blast2table.prl - Parses BLAST reports using BioPerl's Bio::Tools::Blast.pm; writes the data from each HSP in tabular form in a variety of formats.
>> bp_embl2picture.prl - Renders a GenBank or EMBL file into a PNG or GIF image.
>> compare_library.prl - Accepts two files (i and j) containing multiple DNA sequences in FASTA format and compares each sequence in file i to that in file j using a local BLAST installation.
>> count_types.sh - Counts how many files there are of each type in a directory.
>> NCBI_accession_retrieval.sh - Downloads sequences from NCBI in FASTA format when provided with a file containing accession numbers.
>> NCBI_condense_names.prl - Replaces entry names in downloaded FASTA sequences from NCBI with simplier names.
>> NCBI_retrieval.prl - Uses NCBI's Entrez Programming Utilities to perform interactive batch requests to NCBI Entrez.
>> split_fasta.prl - Accepts a file consisting of multiple FASTA formatted sequence records and splits them into multiple files.
>> nanorc.txt - Customized configuration file for use with the GNU Nano 2.0.7 text editor. Allows nucleotide highlighting in FASTA and NEXUS files. Save to your home directory as .nanorc and it will be sourced by Nano at start-up.

FaBox (1.41) - an online fasta sequence toolbox

1. FaBox,
http://onlinelibrary.wiley.com/doi/10.1111/j.1471-8286.2007.01821.x/abstract

FaBox is a collection of simple and intuitive web services that enable biologists to quickly perform typical task with sequence data. The services makes it easy to extract, edit, and replace sequence headers and join or divide data sets based on header information. Other services include collapsing a set of sequences into haplotypes and automated formatting of input files for a number of population genetics and phylogenetic programs, such as ArlequinTCS and MrBayes. The toolbox is expected to grow on the basis of requests for particular services and converters in the future.

2. download
http://users-birc.au.dk/biopv/php/fabox/faq.php

Kent source - bioinformatic operation on fasta and more

1. http://www.biostars.org/p/1852/

2. compile Kent source in Ubuntu
http://genomewiki.ucsc.edu/index.php/Source_tree_compilation_on_Debian/Ubuntu

2013年3月28日星期四

labs in evolutionary genetics

1. barkerlab lab, on genome duplication

http://barkerlab.net/

2Yaniv Brandvain

Population genetics of speciation and mating system evolution


http://yanivbrandvain.wordpress.com/publications/

3. Coop Lab
Population and evolutionary genetics

http://gcbias.org/publications/

visualization of biological data in Google Earth using R2G2, an R CRAN package



Arrigo, N., Albert, L. P., Mickelson, P. G. and Barker, M. S. (2012), Quantitative visualization of biological data in Google Earth using R2G2, an R CRAN package. Molecular Ecology Resources. doi: 10.1111/1755-0998.12012


Set Environment variables - PATH, CLASSPATH, JAVA_HOME, ANT_HOME in Ubuntu


Setting Environment variables in Ubuntu can be tricky. It comes with OpenOffice which requires OpenJDK, so the path for java is set to that of OpenJDK. The command

$ java -version

works but uses the OpenJDK. Once you have installed Sun JDK the location of the JDK should be "/usr/lib/jvm/java-6-sun-1.6.0.20" depending upon the version you have installed. Once you have installed it you may want to update the PATH, CLASSPATH variables and create environment variables like JAVA_HOME, ANT_HOME, etc.

The location to set up the environment variables in Ubuntu is /home/Your_User_Name/.bashrc. You will need to make entries like:

JAVA_HOME=/usr/lib/jvm/java-6-sun-1.6.0.20
ANT_HOME=/home/harkiran/javaTools/apache-ant-1.8.1
PATH=$PATH:$ORACLE_HOME/bin:$JAVA_HOME/bin:$ANT_HOME/bin
CLASSPATH=.:/usr/lib/jvm/java-6-sun-1.6.0.20/lib
export JAVA_HOME
export ANT_HOME
export CLASSPATH
export PATH

The separator to use in Linux between PATH is ":" (colon). Windows uses ";" (semi colon).

Once you have set this up you will also need to create symbolic links in the "/etc/alternatives" directory. You will need administrative privileges to do so.

$ ln -s /usr/lib/jvm/java-6-sun-1.6.0.20/bin/java /etc/alternatives/java
$ ln -s /usr/lib/jvm/java-6-sun-1.6.0.20/bin/javac /etc/alternatives/javac

Once done you can check using the commands java -version, javac -version, ant to check all is working and has been set up properly.
##########################################################
# what I used
# find / -name libjava.so 2> /dev/null
export JAVA_HOME=/usr/lib/jvm/java-6-openjdk-amd64/
export ANT_HOME=/usr/share/ant
export PATH=$PATH:${JAVA_HOME}/bin:${ANT_HOME}/bin
export CLASSPATH=.:/usr/lib/jvm/java-6-sun-1.6.0.20/lib
########
# sudo needed
$ ln -s /usr/lib/jvm/java-6-openjdk-amd64/bin/java /etc/alternatives/java $ ln -s /usr/lib/jvm/java-6-openjdk-amd64/bin/javac /etc/alternatives/javac


2013年3月26日星期二

a bash script that loops through chromosomes

#to write a bash script that loops through chromosomes
#!/bin/bash
for i in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
do
   cbatch "plink --bfile /home/GWAS/GeneralRelease/Imputed/Release3/CleanPlink/GRclean${i} --extract  c${i}.keep1 --maf .05 --out GRclean2 --make-bed"
done