Tuesday, April 30, 2013

Blastx uniref

I chose uniref90 DB from ->
ftp://ftp.ebi.ac.uk/pub/databases/uniprot/uniref/uniref90/
the same can also be found in:
ftp://ftp.uniprot.org/pub/databases/uniprot/uniref/uniref90/

Documentation 

it has 14 205 227 sequences updated on 03/04/2013, Fasta file   6.0 GB

Using blast+
see manual
excellent blast to blast+  

Make blast DB

CMD="dustmasker -in $fasta -infmt fasta -parse_seqids -outfmt maskinfo_asn1_bin -out $fasta_dust.asnb"
$CMD;
echo "$CMD";
(not strictly necesary: identifies and masks low complexity regions of protein, takes at least 4 hours using 10GB)


CMD2="makeblastdb -in $fasta -dbtype prot -parse_seqids -mask_data $fasta_dust.asnb -out uniref90DB_filtered"



blastx -query $FOLD/$filename -db $DB -out $OUTFOLD/$out -outfmt 7 -evalue 1e-10 -num_threads 4

#required resources 8 CPUs 10GB each for 24 hours





Compress / Uncompress

compress tar.gz
tar -pczf try.tar.gz try1/

uncompress tar.gz

 tar -zxvf try.tar.gz

Compress gz

gzip file

Uncompress gz

gunzip -c file.gz > newfile

OR

gzip -d  file.gz > newfile

Select nucleotides from a fasta file using bedtools

module load bedtools 

bedtools getfasta -fi <file.fasta> -bed <file.gff> -fo <outputfile.fasta>

index all fasta files in one line #shell

module load samtools
 
for f in fixed_*fasta; do    samtools faidx $f ; done

How to count nucleotide sequences from a fasta file #python

To count nucleotides in a Fasta file:

from Bio import SeqIO

for seq in SeqIO.parse('chrA01.fasta', 'fasta'):
    s = str(seq.seq)
    print s.count('T') + s.count('G') + s.count('C') + s.count('A')

-------------------------------------------------------
from Bio import SeqIO

for seq in SeqIO.parse('chrA01.fasta', 'fasta'):
    print s.count('T') + s.count('G') + s.count('C') + s.count('A') + s.count('N')

-------------------------------------------------------
from Bio import SeqIO

for seq in SeqIO.parse('chrA01.fasta', 'fasta'):
    print len(seq.seq)