try in a loop #R

If you know there could be and error use try

if it is in a loop and you need to keep going, next!

m <- try(read.table(f, head=FALSE, sep="\t"))
if(class(m) == "try-error") next

compare two files! #shell

Excellent!



#in file 1 but not in file 2 
diff --changed-group-format='%<' --unchanged-group-format='' file1 file2 > output 

#in file 2 but not in file 1
 
diff --changed-group-format='%>' --unchanged-group-format='' file1 file2 > output
 

#Common in both files
diff --changed-group-format='' --unchanged-group-format='%=' file1 file2 > output


From -> linux questions

bedtools

https://github.com/arq5x/bedtools/blob/master/docs/content/general-usage.rst

substring in list #python

if any("abc" in s for s in some_list):
    # whatever

to get all the items containing "abc", use

matching = [s for s in some_list if "abc" in s]

http://stackoverflow.com/questions/4843158/check-in-python-if-a-list-item-contains-a-string 

extract a group of fasta sequences from a file

extract a group of fasta sequences from a file
in python
http://www.biostars.org/p/1709/

from Bio import SeqIO

fasta_file = "fasta_file.fasta" # Input fasta file
wanted_file = "wanted_file.txt" # Input interesting sequence IDs, one per line
result_file = "result_file.fasta" # Output fasta file

wanted = set()
with open(wanted_file) as f:
    for line in f:
        line = line.strip()
        if line != "":
            wanted.add(line)

fasta_sequences = SeqIO.parse(open(fasta_file),'fasta')
with open(result_file, "w") as f:
    for seq in fasta_sequences:
        if seq.id in wanted:
            SeqIO.write([seq], f, "fasta")

grep OR AND etc

OR

grep -E 'pattern1|pattern2' filename

AND

grep -E 'pattern1.*pattern2' filename 

http://www.thegeekstuff.com/2011/10/grep-or-and-not-operators/ 

UniProtKB/Swiss-Prot

 UniprotKB Protein knowledgebase, consists of two sections:

• Swiss-Prot, which is manually annotated and reviewed.
• TrEMBL, which is automatically annotated and is not reviewed.

Includes complete and reference proteome sets.


UniProtKB/Swiss-Prot is the manually annotated and reviewed section of the UniProt Knowledgebase (UniProtKB).
documentation

Downloads > http://www.uniprot.org/downloads

 UniProt Knowledgebase (UniProtKB) Release 
2013_05 (01-May-2013)

UniProtKB/Swiss-Prot: 540,052 entries 

release notes 1st of May 2013 >

 ftp://ftp.uniprot.org/pub/databases/uniprot/relnotes.txt 

http://web.expasy.org/docs/relnotes/relstat.html  

unmodified files, untouched #shell

find ./ -type f -atime +90

grep in loop # shell

In one line

for f in *.excov; do    echo "file $f ";    grep -c 'false' $f; done > lostgenes.txt

split one column into two by a string #R

library(stringr)
str_split_fixed(column, "string", 2)

call R files from shell

module load R

R --no-save < file.R

Sum matrices from a list #R

Mergefiles <- function(FOL, listM){
        path <- getwd()
        setwd(FOL)
        print(length(listM))
        newmatrix <- Reduce('+', listM)
        #add first column names
        newmatrix[,1]<- names(newmatrix)[-1]
        #write file
        write.table(newmatrix,file="mergedmatrixXA.txt",sep="\t", row.names=FALSE, quote=FALSE)
        setwd(path)
}

UniProt Go

UniProt GOA (Uniprot GO Annotation)


http://www.ebi.ac.uk/GOA
Readme file

UniProt (Swiss-Prot and TrEMBL)

gene_association.goa_uniprot 

Locations: 
ftp://ftp.geneontology.org/.../gene_association.goa_uniprot.gz ftp://ftp.ebi.ac.uk/.../gene_association.goa_uniprot.gz 

This file contains all GO annotations for proteins in the UniProt KnowledgeBase (UniProtKB).

Read me file ftp://ftp.ebi.ac.uk/pub/databases/GO/goa/UNIPROT/gp_association_readme 

GAF 2.0 file format

 colNames <- c("DB",  "DB_Object_ID",  "DB_Object_Symbol",  "Qualifier",  "GO_ID",  "DB:Reference",  "Evidence",  "With",  "Aspect",  "DB_Object_Name",  "Synonym",  "DB_Object_Type",  "Taxon_ID",  "Date",  "Assigned_By",  "Annotation_Extension",  "Gene_Product_Form_ID")

17 columns

gene_association.goa_arabidopsis 

Locations:
ftp://ftp.geneontology.org/.../gene_association.goa_arabidopsis.gz ftp://ftp.ebi.ac.uk/.../gene_association.goa_arabidopsis.gz 

This file contains the GO assignments for the proteins of the arabidopsis UniProtKB Complete 
Proteome. 

Read me file ftp://ftp.ebi.ac.uk/pub/databases/GO/goa/ARABIDOPSIS/README 

Extras

might search for individual IDs in the uniref90DB
http://www.ebi.ac.uk/ebisearch/search.ebi?db=uniref90&t=B9DGW5

I might check Arabiposis annotation from TAIR10
ftp://ftp.arabidopsis.org/home/tair/Ontologies/Gene_Ontology/ATH_GO_GOSLIM.txt

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)