text:
How can I read FCS files using open source libraries?
FCS is a patented data format used for storing flow cytometry data. The most recent version is FCS3.1. There is some documentation on the format, but there is no information on how to read these files. There are some R packages and a MATLAB code to read an FCS file, but I am looking for standard libraries developed either by the FCS consortium or any other group. I also wish to know if FCS is a subset of an existing standard data format that can be read by a standard library using any programming language.
Finally, I would want to convert these files to an easily readable format like HDF5.

I find strange you ask for standard libraries but you exclude the existing standard tools.

@Llopis the OP isn't rejecting standard tools, on the contrary, he is looking for standard libraries.

@terdon I am not sure how standard libraries are decided but a R/Bioconductor package or a MATLAB one seems to me quite standard (Except if the FCS consortium provides one of course).

@Llopis By standard, I mean is there an OS/language independent interface that allows reading of these files. More importantly, I want to know how FCS files can be read by any general language. I am expecting something like a pseudocode, for this purpose.

I’m still not really sure what you want: the documentation you link to is as close to a general pseudocode as you will get: it contains a detailed, concise specification of the format. That’s the same information that a pseudocode would convey, and it’s the conventional way in which data formats are described (pseudocode is virtually never used for this). That, in addition to the existence of libraries for common analysis languages (R, Python, Matlab), should be sufficient to work with the format.

@WYSIWYG In addition to what Konrad said, the flowCore in  Bioconductor is available for Windows, iOS and Linux.

@KonradRudolph What I want is an instruction on how these files are to be read. The documentation describes the format but I could not figure out the general procedure to read the format which is what I meant by pseudocode. Just assume that I do not want to use bioconductor;  lets say that I want to read these files in C. How should I go about that? Do the FACS machines use a standard mechanism/io library for reading/writing FCS files or they use their own codes? I assumed that if a format is defined by a consortium, they would have provided a standard interface for rw operation.

Why do you say .fcs is a patented format? The disclaimer in the linked documentation just states "Attention is called to the possibility that implementation of this standard may require use of subject matter covered by patent rights".
It may not even be possible to patent it, see: https://patents.stackexchange.com/questions/11968/can-patents-on-computer-file-formats-be-obtained

A few years ago, I wrote a python script to convert FCS files into tab-separated format. It was far from handling all the possibilities that the format description offers, but at least it worked for some of the files produced on one of our machine: http://www.igh.cnrs.fr/equip/Seitz/en_equipe-programmes.html
The format documentation I found enabled decoding (see section 3 of the pdf you mention), but it requires reading data in binary mode.
The general idea of this format (and, I guess, many other binary formats), is that there is a header zone at the beginning of the file with a defined number of fields encoding numbers indicating how the rest of the file is structured. So a first phase is to parse this header, following the description given in the documentation of the format. The information extracted from the header tells where to find the data and how it is encoded, still according to rules described in the format documentation.
In case this may be useful, and for the record, here is the code from the above-mentioned script (after stripping comments, some of which are merely copied from the format documentation, and adding a few ones):
#!/usr/bin/env python
"""This script tries to read FCS flow cytometry data.
Format parsing inspired by information found here:
http://isac-net.org/Resources-for-Cytometrists/Data-Standards/Data-File-Standards/Flow-Cytometry-Data-File-Format-Standards.aspx
"""

import re
# To decode binary-encoded data
import struct
import sys

class Parameter(object):
    """This object represents one of the parameter types that are present in a DATA segment of a FCS file."""
    __slots__ = ("p_name", "p_bits", "p_range", "p_ampl", "parser")
    def __init__(self, p_name, p_bits, p_range, p_ampl):
        self.p_name = p_name
        self.p_bits = p_bits
        self.p_range = p_range
        self.p_ampl = p_ampl
        # Function for parsing a value of the parameter in the data segment
        self.parser = None

##############################################
# Here starts the parsing of the header part #
# which tells where the other parts are.     #
##############################################

f = open(sys.argv[1], "rb")
# The format name is encoded in 6 letters
# An ASCII letter is coded with one octet
file_format = "".join([f.read(1) for __ in range(6)])
sys.stdout.write("Format: %s\n" % file_format)
# The format descriptions reserves 4 octets that we skip
skip = f.read(4)
# 8 octet chunks encode the start and end positions
# of different parts of the data
text_start = int(f.read(8).strip(" "))
text_end = int(f.read(8).strip(" "))
data_start = int(f.read(8).strip(" "))
data_end = int(f.read(8).strip(" "))
analysis_start = int(f.read(8).strip(" "))
analysis_end = int(f.read(8).strip(" "))

if (analysis_start and analysis_end):
    sys.stderr.write("Cannot deal with ANALYSIS segment of an FCS file.\n")

####################################################
# Here starts the parsing of the "TEXT" portion    #
# which describes how the data proper is organized #
####################################################
f.seek(text_start)
# The first character in the primary TEXT segment is the ASCII delimiter character.
sep = f.read(1)
if sep not in ["_", "@"]:
    alt_sep = "_@_"
elif sep not in ["_", "|"]:
    alt_sep = "_|_"
else:
    assert sep not in ["+", "|"]
    alt_sep = "+|+"
text_segment = f.read(text_end - text_start)

fields = text_segment.split(sep)

info = {}

i = 0
while i < len(fields) - 1:
    key = fields[i]
    i += 1
    val = fields[i]
    i += 1
    # Keywords are case insensitive, they may be written in a file in lower case, upper case, or a
    # mixture of the two. However, an FCS file reader must ignore keyword case. A keyword value may
    # be in lower case, upper case or a mixture of the two. Keyword values are case sensitive.
    info[key.upper()] = val
print "%s events were detected." % info["$TOT"]
print "Each event is characterized by %s parameters" % info["$PAR"]

if info["$NEXTDATA"] != "0":
    sys.stderr.write("Some other data exist in the file but hasn't been parsed.\n")

# L - List mode. For each event, the value of each parameter is stored in the order in which the
# parameters are described. The number of bits reserved for parameter 1 is described using the
# $P1B keyword. There can be only one set of list mode data per data set. The $DATATYPE
# keyword describes the data format. This is the most versatile mode for the storage of flow
# cytometry data because mode C and mode U data can be created from mode L data.
assert info["$MODE"] == "L"

parameters = []

# indices of the parameters
p_indices = range(1, int(info["$PAR"]) + 1)
for i in p_indices:
    p_name = info["$P%dN" % i]
    p_bits =  info["$P%dB" % i]
    p_range = info["$P%dR" % i]
    p_ampl =  info["$P%dE" % i]
    parameters.append(Parameter(p_name, p_bits, p_range, p_ampl))

sys.stdout.write("The parameters are:\n%s\n" % "\t".join([par.p_name for par in parameters]))

# How are 32 bit words organized
if info["$BYTEORD"] == "4,3,2,1":
    endianness = ">"
else:
    endianness = "<"
    assert info["$BYTEORD"] == "1,2,3,4"

    # I stripped a long comment which is just a copy of the documentation
# Type of data:
if info["$DATATYPE"] == "I":
    for par in parameters:
        nb_bits = int(par.p_bits)
        assert nb_bits % 8 == 0
        nb_bytes = nb_bits / 8
        # Determine format string for unpacking (see https://docs.python.org/2/library/struct.html)
        if nb_bytes == 1:
            c_type = "B" # unsigned char
        elif nb_bytes == 2:
            c_type = "H" # unsigned short
        elif nb_bytes == 4:
            c_type = "L" # unsigned long
        elif nb_bytes == 8:
            c_type = "Q" # unsigned long long
        else:
            raise ValueError, "Number of bytes (%d) not valid for an integer (see https://docs.python.org/2/library/struct.html#byte-order-size-and-alignment)." % nb_bytes
        fmt = "%s%s" % (endianness, c_type)
        p_range = int(par.p_range)
        def parser(data):
            value = struct.unpack(fmt, data.read(nb_bytes))[0]
            try:
                assert value < p_range
            except AssertionError:
                print "Value %s higher than %d" % (str(value), p_range)
            return value
        par.parser = parser
    pass
else:
    raise NotImplementedError, "Only the parsing of integer value has been implemented so far."

out_file = open(sys.argv[2], "w")
out_file.write("#amplification_types\t" + "\t".join([par.p_ampl for par in parameters]) + "\n")
out_file.write("parameters\t" + "\t".join([par.p_name for par in parameters]) + "\n")
i = 1
##############################################
# Here starts the parsing of the data proper #
##############################################
f.seek(data_start)
while f.tell() < data_end:
    values = []
    for par in parameters:
        values.append(par.parser(f))
    out_file.write("%d\t" % i + "\t".join(map(str, values)) + "\n")
    i += 1
out_file.close()
f.close()

The general idea of this format (and, I guess, many other binary formats), is that there is a header zone at the beginning of the file with a defined number of fields encoding numbers indicating how the rest of the file is structured. So I begin by parsing this header, following the description given in the documentation of the format. The information extracted from the header tells me where to find the data and how it is encoded, still according to rules described in the format documentation. I'm not good at writing pseudocode, but I can try to add comments in the python code.

I maintain several Python libraries for interacting with flow cytometry data (FCS files). The most recent of which, FlowKit, provides functionality to export the event data as a NumPy array, Pandas DataFrame, or CSV text. Using the h5py library, it's trivial to export as an HDF5 file

R/Bioconductor has a number of different flow cytometry processing packages. One place to start for looking at cytometry data from a high level would be openCyto (or its vignette), which is a large set of tools for basic extraction and analysis of FCS files.
I have looked in the past at the FCS files as an R structure using flowCore. Loading a single FCS file is fairly straightforward and follows a familiar R pattern:
file.name <- "/dir/file.fcs"
x <- read.FCS(file.name, transformation=FALSE)
summary(x)

Asking about converting FCS files into an "easily readable format like HDF5" doesn't seem like the right question. HDF5 is a container format and shares a lot of similarity to file systems. I've found it best to keep FCS files as they are, as it is a compact, standardised binary format.

HDF5 is not a major concern (I just like it). I am still curious about how something like read.FCS works. I believe if FCS is a standard format then there should be standard (non-proprietary) libraries to read/write FCS files. As of now I'm forced to use R or Matlab (and still there are no standard libraries here too. I guess there are just different packages/programs made by different groups/individuals; no single standard). What if I want to analyse a FCS file on C, Fortran or Python (etc)?


text:
A tool or webserver for building PSSM matrix
I have some protein sequences and I want to build a position-specific scoring matrix (PSSM) for them and then upload this PSSM to NCBI PSI-BLAST. I used CHAPS program for this pupose but uploading the output PSSM gave me an error in NCBI PSI-BLAST. Do you know any tool or webserver for getting PSSMs for a group of protein sequences which then can work in NCBI PSI-BLAST?
This is a part of my PSSM file:
PssmWithParameters ::= {
pssm {
    isProtein TRUE ,
    numRows 28 ,
    numColumns 131 ,
    byRow FALSE ,
    query
    seq {
        id {
            other {
                accession "WP_000208753" } } ,
        inst {
            repr raw ,
            mol aa ,
            length 131 ,
            seq-data
            ncbieaa "MTTKRKPYVRPMTSTWWKKLPFYRFYMLREGTAVPAVWFSIELIFGLFALKNGPEAW
AGFIDFLQNPVIVIINLITLAAALLHTKTWFELAPKAANIIVKDEKMGPEPIIKSLWAVTVVATIVILFVALYW" } } ,
    intermediateData {
        freqRatios {
            { 0, 10, 0 } ,
            { 564418841, 10, -10 } ,
            { 0, 10, 0 } ,
            { 11768571, 10, -9 } ,
            { 185838265, 10, -10 } ,
            { 31496547, 10, -9 } ,
            { 3872857, 10, -8 } ,
            { 291750464, 10, -10 } ,
            { 128450763, 10, -10 } ,
            { 759856221, 10, -10 } ,
            { 359173937, 10, -10 } ,
            { 179865517, 10, -9 } ,
            { 14537895, 10, -8 } ,
            { 212921456, 10, -10 } ,
            { 220554141, 10, -10 } ,
            { 368516324, 10, -10 } ,
            { 319343525, 10, -10 } ,
            { 426173953, 10, -10 } ,
            { 463659523, 10, -10 } ,
            { 817322186, 10, -10 } ,
            { 811693964, 10, -11 } ,
            { 0, 10, 0 } ,
            { 227810064, 10, -10 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 768325726, 10, -10 } ,
            { 0, 10, 0 } ,
            { 1425562, 10, -8 } ,
            { 372685285, 10, -10 } ,
            { 466727421, 10, -10 } ,
            { 185741084, 10, -10 } ,
            { 427328646, 10, -10 } ,
            { 122594965, 10, -10 } ,

Which error did it give you the NCBI PSI-BLAST? Could you paste here the image with the options used? What have you tried to overcome this problem?

The error was: Message ID#9 Error: PSSM processing error: Incomplete PSSM with Parameters provided on input: missing score matrix

So it seems that you need to provide the PSSM in certain format and CAPS didn't manage to create it in the expected format. Could you post the PSSM here (or in an accessible way) ? Now that I look in the help I couldn't find any reference to a specific format (just that it should be from their own output ...)

I do not exactly know how to post PSSM file here. Could you please guide me how to load a file here?

You can't but you can copy paste the content or add a link to dropbox or a similar server

@Llopis it should be fine posted here as code.

@Sara see formatting tools for help on formatting your posts. For things like this, you want to use code ({}) which uses a fixed-width font.

I think you could use standalone blast, one example here https://www.biostars.org/p/15482/

If you take a look at my answer in this BioStars post, you can generate a PSSM using AlignIO in Biopython:

You can use biopython's code for generation of PSSM.
Here is the link for the source code which can help you.

I'm not 100% sure what you're trying to do with regards to your search in the first place, but I recommend using Sean Eddy's program Jackhammer it's similar to PSI-BLAST but uses Hammer HMM underpinnings.
As input it can take a single sequence, a multiple sequence alignment, or an existing HMM and then allow you to search in an iterative model building way.  You can also download an HMM too at any stage in the iterative search, as well as the results / alignment.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Obtaining data table headers from GEO using GEOquery
For a study in GEO, I would like to obtain the data table header descriptions, specifically the "VALUE" column for all samples in the study.
If you go here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE99511)
and then scroll down and click one of the samples: let's choose "GSM2644971". Then scroll down and you should see "Data table header descriptions" and below that you should see "VALUE Normalized (provided the normalization method) Average Beta". That information is what I want.
I tried using assayData() from the Biobase package but I didn't know whether the method takes a sample, a matrix of samples, or something else as a parameter. Apparently it takes a S4 and then returns some random text but I don't know how to use that random text.
EDITED:
For example, if I do:
library(Biobase)
library(GEOquery)

getgeo<-getGEO("GSE99511")
assayData(getgeo$GSE99511_series_matrix.txt.gz)

When I use assayData() I would want it to return "Normalized (provided the normalization method) Average Beta".
Instead it returns:
<environment: 0x110674178>

Please let me know if the question doesn't make sense or if there's another method that I should be using.

It would be helpful if you could add a specific example we can copy and test on our own machines. Show us some simple data, or a way of generating appropriate random data and then show the output you get and the output you want to get instead.

you can download the GEO metadata directly from the GEO page...look for the SOFT file or the series matrix files, for example.

The assayData requires a SummarizedExperiment, if you want to download the description of the value you need to look at the GEOQuery package when I can (probably in 24h )

I want to download it through GEOquery so that I can get data faster. I also updated my question to add more detail.

Did you check what return getGEO, and what do you need to extract the data you want?

https://stackoverflow.com/questions/44856211/retrieving-data-table-headers-from-geo-using-geoquery

The short answer is that if you are seeing information in a GEO sample, then it is the sample that you need to download and access using the appropriate method. In this case, that method is Columns():
library(GEOquery)
gsm <- getGEO("GSM2644971")
Columns(gsm)

          Column      Description
1         ID_REF
2         VALUE       Normalized (provided the normalization method) Average Beta
3         Detection   Pval

The longer answer: you will really benefit from sitting down with the GEOquery documentation and trying to understand the objects that it uses and their methods. It is not something you can work out by trial and error for every case-by-case basis. The tutorial is a really good starting point (you will find Columns() in there, for example).

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
PASA pipeline: compare experimental transcripts to the reference annotation
I would like to ask if anyone has experience in running a subset of the PASA pipeline, in particular for the reconciliation of some experimental 'transcripts' with the reference annotation.
In more detail, I am working with RNA-seq data from D. melanogaster. I have reconstructed the 'transcripts' using Trinity. I have aligned these 'transcripts' to the reference genome using GMAP. Now I would like to match these 'experimental transcripts' with the reference annotation, to see how they compare.
I was wondering if I can run just a subset of the whole PASA pipeline, so basically skipping the step of alignment to the reference and using the alignment file that I generated externally.
Also, I would be glad to receive more general redirecting to any released software allowing the comparison of experimental transcripts to the reference annotation.
Cross-posted on the PASA Google Group pasapipeline-users
Cross-posted on biostars

Go easy on the cross-posting. You've already directly asked the authors, just be patient for their reply.

Hi Amina Echhiki, thanks for your question and welcome to Bioinformatics Stack Exchange. As you've included this (or a similar question) on other sites, it would be appreciated if you could report on the answers you have received as your own answer to this question. Stack Exchange encourages people to answer their own questions as they gain knowledge about the subject.

Answer from the authors:

PASA ends up using gff3 format instead of sam format for uploading the
  alignments.  If somebody has the gmap gff3 output, then can be done to
  just upload that directly using the custom alignment importer, but not
  with sam. They suggest just having PASA rerun GMAP as part of its
  regular routine.


text:
Workflow of metabarcoding analysis
For a diet analysis of an insect-eating animal, all species in a sample shall be identified. For this a sequencing of the metagenomic sample was done, where the COI/COX region was used as a barcode and amplified. This resulted in 20 million reads (after quality control).
While looking for an apropriate workflow, I am faced with these issues:

For many tools 20 million reads are to much to handle (BLAST for
instance)
For the majority of metagenomic tools the focus is on micro-organisms, not insects
The database to search against must contain COI/COX data, mainly from insects, however, the animal could also eat snails or frog eggs...
Since the reads are amplicons some analysis are not recommended

Leaving me with these questions:

What are the typical analysis steps in a workflow like this?
Which tools could I use for my purpose?
Where could I get a database of COI/COX data?

I'm working on something similar as a side project. I'll detail the steps that I'm using, but I'm curious if others can suggest a better way.

You don't mention if your data are paired-end, but in my case I have paired-end sequences that overlap in the middle. I used FLASh to merge the reads and seqtk to convert the reads to fasta format.
I used the cluster_fast algorithm in vsearch to cluster sequences at 100% identity, then filtered out rare sequence types (<10 sequences) to account for sequencing errors.
This resulted in a manageable number of unique sequences, which I searched against a custom database using vsearch, but which could also have been searched against the nr database using blast.

Regarding your third question ("Where could I get a database of COI/COX data?"), you can apparently search and download fasta sequences from BOLD.
However, I don't see a "select all" button to select all records found by one search in a single fasta file.
There seems to exist an R API: https://cran.r-project.org/web/packages/bold/bold.pdf which refers to this: http://www.boldsystems.org/index.php/resources/api#sequenceParameters>.
From this I seem to understand that one can directly query using urls:
wget -O Gastropoda.fa http://www.boldsystems.org/index.php/API_Public/sequence?taxon=Gastropoda&marker=COI

And I get a fasta file containing COI sequences for Gastropoda. There may be ways to refine the query geographically, using the geo parameter.

Yes, I already downloaded some of the BOLD data. I circumvented the missing "select all" button by using a firefox addon for selection boxes. I did not know that there is an R API, that will be helpful.


text:
Obtaining all protein sequences with a particular domain architecture from Pfam
I want to get the alignment of chain A of 1kf6 (PDB ID) from the pfam database here. This protein chain has two main domains (FAD_binding_2 and    Succ_DH_flav_C). In pfam there is a link to one of these domains and after clicking in one of the mentioned domains in the table below the page, another page comes at the top of which there is a link to architectures.
For example, if we click on "FAD_binding_2" in the table, there are 220 architectures that have this domain and after clicking on that button (at the top of page), a page opens here. In this page 7471 sequences with the following architecture: FAD_binding_2, Succ_DH_flav_C, resembles most to my chain (because it has both FAD_binding_2 and Succ_DH_flav_C domains).
If we want to see all 7471 sequences with that architecture we can click on the show button under that architecture and finally we can see all of that 7471 architecture resembling our chain A sequence. Down the page we can see other architectures that all have FAD_binding_2, but some of them do not have both FAD_binding_2 and Succ_DH_flav_C domains (or some of them have other extra domains that we are not interested in) and so we do not need their sequences and their alignments.
When we want to get the alignment of full sequences at the left side of the page, we can see that we can just get the alignment of all 15696 sequences and it seems that there is no way to get the alignment of our favourite 7471 sequences. I would like to know if there is any way to get the alignment or at least the fasta file of our favourite 7471 sequences (which have both of our favourite domains not just one of them)?

I believe you can do this using SQL files of the database.
First you need to get the integer accession of this architecture. architecture table has the following columns:

auto_architecture (integer id, this is what we need)
architecture (a string describing the architecture)
type_example (example sequence)
no_seqs (number of sequences)
architecture_acc (accessions of the architectures)

You need the following line:
3900719357      FAD_binding_2~Succ_DH_flav_C    Z9JRB3  7471    PF00890 PF02910

Now you need to get the sequences with this architecture id. First, download the file with all the sequences (warning, this one is more than 7Gb). With this one you need column #1 (pfam sequence accesion), #12 (sequence itself) and #16 (architecture integer id).
If you do not want to import files in the database, you can use awk to export all the sequences to fasta format.
zcat pfamseq.txt.gz | awk -F\\t '{ if ($16=="3900719357") print ">" $1 "\n" $12}' > sequences.fasta

Of couse, if you need to perform multiple queries like this, it is much easier to import tables into the database. Then you can perform both queries in a single SQL operation (assuming you know the architecture string):
SELECT pfamseq_acc, sequence
    FROM architecture
    JOIN pfamseq ON pfamseq.auto_architecture=architecture.auto_architecture
    WHERE architecture="FAD_binding_2~Succ_DH_flav_C";

(I haven't tried this SQL query, so maybe you need to correct it a bit)

I think the idea of not importing 7Gb database is better because I have a lot of protein sequences that I want to do similar task on. I prefer trying the second idea, but I don't know why I couldn't get the integer accession of this architecture from above "architecture table" link

Not sure I understand your comment well. In my example I assume you only have the architecture string. If you already have your auto_architecture accession, you can simply write SELECT pfamseq_acc, sequence FROM pfamseq WHERE auto_architecture=3900719357.

An alternative way is using the Hmmer website, and do a hmmsearch with the Pfam accession, forex. (PF02910). Navigate to the domain tab and find the domain architecture that you are interested in and press view scores. That way, only the sequences you are interested in will be downloadable under the tab downloads.
Step-by-Step guide:
1- Go to http://www.ebi.ac.uk/Tools/hmmer/
2- Search with hmmsearch http://www.ebi.ac.uk/Tools/hmmer/search/hmmsearch
3- Use Accession search and type in PF02910
4- Navigate to Domain
5- To the right side of
9742 sequences with domain architecture: FAD_binding_2, Succ_DH_flav_C
Press View scores.
6- It should show (Your results have been filtered) and Query matches (9742)
7- Navigate to Download and download the sequences in the format of interest.
8- Repeat for domain architectures of interest.
I hope you find this helpful.
Kind Regards,
Sara

That does not appear to be possible using PFAM. I would suggest trying a different database. I managed to perform the same query using SMART (http://smart.embl-heidelberg.de/), via entering "Pfam:FAD_binding_2 AND Pfam:Succ_DH_flav_C" under "Domain selection" in the "Architecture Analysis" section of the search screen, which resulted in 9921 protein sequences with those two domains as you can see here.
On the results screen you can select all the sequences and under "Action" at the top of the screen choose "Download Protein Sequences as FASTA file".

It's a good idea but the problem is that it's not specific. I need the proteins that just have these two domains and no more. In these 9921 protein sequences, there are sequences which have other domains that make my search inspecific

you gave me a good idea. although it's inspecific but I could simply make it specific. First in "action selection" part, I chose "download protein sequences as a fasta file" and downloaded the fasta file of all sequences and then I chose "generate a Newick tree and iTol formatted dataset" . Then in "iTOL protein dataset file" part, I clicked on "download as tab delimited plain text" and downloaded all protein identifiers with the name of their pfam domain

Then I used python for deleting protein identifiers which had other domains. After that I could go to the fasta file and delete the fasta sequence of extra proteins


text:
Selecting sites from VCF which have an alt AD > 10
I have high-depth variant calling created using the HaplotypeCaller with --output_mode EMIT_ALL_SITES I'm interested in finding all sites (regardless of genotype call heterozygous or homozygous) where at least one of the alternative alleles have an AD value (Allelic Depth) greater than 10, I.e. are supported by more than 10 reads. Also ideally I want back more than just the first alternative allele.  Note that I don't want back lines of VCF were we only see an AD count for the ref allele only.
So in the example VCF snippet below I'm wanting to select lines: 6,7,8,12,13 and 14, which have GT:AD values 1/1:1,988:989 0/1:116,92 0/1:220,234 0/1:62,611 1/1:0,109 respectively.
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  12908_DIAG
3   187446740   .   T   .   Infinity    .   AN=2;DP=1095;MQ=60.00   GT:AD:DP    0/0:1095:1095
3   187446741   .   C   .   Infinity    .   AN=2;DP=1117;MQ=60.00   GT:AD:DP    0/0:1117:1117
3   187446752   .   A   .   Infinity    .   AN=2;DP=1297;MQ=60.00   GT:AD:DP    0/0:1297:1297
3   187446763   .   C   .   Infinity    .   AN=2;DP=1494;MQ=60.00   GT:AD:DP    0/0:1494:1494
3   187451574   .   C   .   Infinity    .   AN=2;DP=1493;MQ=60.00   GT:AD:DP    0/0:1493:1493
3   187451609   rs1880101   A   G   39794.03    .   AC=2;AF=1.00;AN=2;BaseQRankSum=1.859;ClippingRankSum=0.000;DB;DP=995;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;MQRankSum=0.000;QD=24.56;ReadPosRankSum=0.406;SOR=8.234  GT:AD:DP:GQ:PL  1/1:1,988:989:99:39808,2949,0
4   1803279 .   T   G   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-6.652;ClippingRankSum=0.000;DP=245;ExcessHet=3.0103;FS=89.753;MLEAC=0;MLEAF=0.00;MQ=59.97;MQRankSum=0.000;ReadPosRankSum=-2.523;SOR=6.357   GT:AD:DP:GQ:PL  0/0:211,23:234:99:0,364,6739
4   1803307 rs2305183   T   C   2486.60 .   AC=1;AF=0.500;AN=2;BaseQRankSum=-5.049;ClippingRankSum=0.000;DB;DP=215;ExcessHet=3.0103;FS=1.110;MLEAC=1;MLEAF=0.500;MQ=59.97;MQRankSum=0.000;QD=11.95;ReadPosRankSum=-0.045;SOR=0.809  GT:AD:DP:GQ:PL  0/1:116,92:208:99:2494,0,3673
4   1803671 .   C   A   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-0.880;ClippingRankSum=0.000;DP=450;ExcessHet=3.0103;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;MQRankSum=0.000;ReadPosRankSum=-0.953;SOR=0.572    GT:AD:DP:GQ:PL  0/0:445,2:447:99:0,1272,15958
4   1803681 .   T   C   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-1.654;ClippingRankSum=0.000;DP=483;ExcessHet=3.0103;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;MQRankSum=0.000;ReadPosRankSum=-0.422;SOR=0.664    GT:AD:DP:GQ:PL  0/0:479,2:481:99:0,1408,18538
4   1803703 .   A   G   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-1.704;ClippingRankSum=0.000;DP=458;ExcessHet=3.0103;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;MQRankSum=0.000;ReadPosRankSum=0.299;SOR=0.497 GT:AD:DP:GQ:PL  0/0:454,2:456:99:0,1325,18095
4   1803704 rs2234909   T   C   6676.60 .   AC=1;AF=0.500;AN=2;BaseQRankSum=-2.605;ClippingRankSum=0.000;DB;DP=456;ExcessHet=3.0103;FS=1.753;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=14.71;ReadPosRankSum=0.324;SOR=0.849   GT:AD:DP:GQ:PL  0/1:220,234:454:99:6684,0,6366
4   1803824 rs2305184   C   G   2030.60 .   AC=1;AF=0.500;AN=2;BaseQRankSum=8.083;ClippingRankSum=0.000;DB;DP=124;ExcessHet=3.0103;FS=6.128;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=16.51;ReadPosRankSum=0.180;SOR=0.096    GT:AD:DP:GQ:PL  0/1:62,61:123:99:2038,0,1766
4   1805296 rs3135883   G   A   3876.03 .   AC=2;AF=1.00;AN=2;DB;DP=110;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=29.22;SOR=9.401    GT:AD:DP:GQ:PL  1/1:0,109:109:99:3890,326,0

A dropbox link for file
I'd initially considered using GATK's SelectVariants but I'm not sure JEXL has the ability to select out what I want specifically other than a blanket AD > 10 which will give me both ref and alt alleles with AD > 10.  Perhaps there is a bioawk solution or something more elaborate with coreutils which could successfully return sites with an alt AD count > 10?

cross-posted: http://gatkforums.broadinstitute.org/gatk/discussion/comment/40003#Comment_40003

This now works with the development version of Bcftools v1.5 (commit 4f134df). Thanks to Petr Danecek for adding the feature. I expect this feature to make its way into the next release of Bcftools:
git clone git://github.com/samtools/htslib.git
git clone git://github.com/samtools/bcftools.git
(cd bcftools; make)

bgzip Test.vcf
./bcftools/bcftools index Test.vcf.gz
./bcftools/bcftools filter -i 'AD[0:1-] > 10' Test.vcf.gz

Output without header (I have modified the second line to be tri-allelic to demonstrate the filtering works):
3   187451609   rs1880101   A   G   39794   PASS    AC=2;AF=1;AN=2;BaseQRankSum=1.859;ClippingRankSum=0;DB;DP=995;ExcessHet=3.0103;FS=0;MLEAC=2;MLEAF=1;MQ=60;MQRankSum=0;QD=24.56;ReadPosRankSum=0.406;SOR=8.234   GT:AD:DP:GQ:PL  1/1:1,988:989:99:39808,2949,0
4   1803279 .   T   G   0   PASS    AC=0;AF=0;AN=2;BaseQRankSum=-6.652;ClippingRankSum=0;DP=245;ExcessHet=3.0103;FS=89.753;MLEAC=0;MLEAF=0;MQ=59.97;MQRankSum=0;ReadPosRankSum=-2.523;SOR=6.357 GT:AD:DP:GQ:PL  0/0/0:211,3,34:234:99:0,364,6739
4   1803307 rs2305183   T   C   2486.6  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-5.049;ClippingRankSum=0;DB;DP=215;ExcessHet=3.0103;FS=1.11;MLEAC=1;MLEAF=0.5;MQ=59.97;MQRankSum=0;QD=11.95;ReadPosRankSum=-0.045;SOR=0.809   GT:AD:DP:GQ:PL  0/1:116,92:208:99:2494,0,3673
4   1803704 rs2234909   T   C   6676.6  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-2.605;ClippingRankSum=0;DB;DP=456;ExcessHet=3.0103;FS=1.753;MLEAC=1;MLEAF=0.5;MQ=60;MQRankSum=0;QD=14.71;ReadPosRankSum=0.324;SOR=0.849  GT:AD:DP:GQ:PL  0/1:220,234:454:99:6684,0,6366
4   1803824 rs2305184   C   G   2030.6  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=8.083;ClippingRankSum=0;DB;DP=124;ExcessHet=3.0103;FS=6.128;MLEAC=1;MLEAF=0.5;MQ=60;MQRankSum=0;QD=16.51;ReadPosRankSum=0.18;SOR=0.096    GT:AD:DP:GQ:PL  0/1:62,61:123:99:2038,0,1766
4   1805296 rs3135883   G   A   3876.03 PASS    AC=2;AF=1;AN=2;DB;DP=110;ExcessHet=3.0103;FS=0;MLEAC=2;MLEAF=1;MQ=60;QD=29.22;SOR=9.401 GT:AD:DP:GQ:PL  1/1:0,109:109:99:3890,326,0

Edited (Dec 2020): bcftools filter syntax has been updated, answer now reflects changes found in version v1.11 (not sure at what point these were introduced).  The old query 'AD[1-] > 10' would produce the error The FORMAT tag AD can have multiple subfields, run as AD[sample:subfield]. A discussion of these changes can be found here.

So just revisiting this and the query now errors in current versions of bcftools: The FORMAT tag AD can have multiple subfields, run as AD[sample:subfield]

Feel free to edit the post. Anyone can...

Also, thanks for the update!

Okay I think the answer is now: 'AD[0:1-] > 10' just testing...

using vcfilterjs
and the following script:
function accept(vc)
    {
    var i,j;
    for(i=0;i< vc.getNSamples();++i)
        {
        var genotype = vc.getGenotype(i);
        if(!genotype.hasAD()) continue;
        var ad = genotype.getAD();
        /* loop over AD starting from '1' =  first ALT */
        for(j=1  ;j< ad.length ;++j)
            {
            if(ad[j]>10) return true;
            }
        }
    return false;
    }
accept(variant);

usage:
java -jar dist/vcffilterjs.jar -f script.js Test.vcf

I think this system you've created looks to be very flexible I'd be very interested in seeing some more walkthroughs and documented examples, on Github because I'm sure it can solve a lot of complex variant filtration issues.

@MatthewBashton thanks ! and there is a much more faster alternative version since today (not javascript , but 100% java based ) http://lindenb.github.io/jvarkit/VcfFilterJdk.html

@DanielStandage not braces, just lambdas using vcffilterjdk : java -jar dist/vcffilterjdk.jar  -e 'return variant.getGenotypes().stream().filter(G->G.hasAD() && java.util.Arrays.stream(G.getAD()).skip(1).filter(AD->AD>10).findAny().isPresent()).findAny().isPresent();' input.vcf

Nice I'll have to give the jdk version a go too.  I have no issues your braces or indentation style ¯\(ツ)/¯

You can do this in Hail:
from hail import *
hc = HailContext()
(hc.import_vcf('test.vcf')
 .filter_variants_expr('gs.exists(g => g.ad[1:].exists(d => d > 10))')
 .export_vcf('filtered.vcf'))

This works with any number of samples and will keep the variants where at least one sample has a genotype with an alternate allele support by more than 10 reads.
To verify we got the expected 6 variants:
>>> hc.import_vcf('filtered.vcf').count()
(1L, 6L)

count returns the number of samples (1) and number of variants (6).
Take a look at the getting started page or tutorials if you want to try it out!

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
How to append numbers only on duplicates sequence names?
I have a reference database with contains 100s of sequences in fasta format. Some of these sequences have duplicate names like so:
>1_uniqueGeneName
atgc
>1_anotherUniqueGeneName
atgc
>1_duplicateName
atgc
>1_duplicateName
atgc

Is is possible to run through a large file like this and change the names of only the duplicates?
>1_uniqueGeneName
atgc
>1_anotherUniqueGeneName
atgc
>1_duplicateName_1
atgc
>1_duplicateName_2
atgc

Do you want to only change the names of the duplicates, I mean of the second, third etc seqs with the same name or also of the 1st as you have done? I mean, do you really need duplicate_1, duplicate_2 or is duplicate, duplicate_1 enough? The latter is far simpler since the former will require you to read the file twice.

Cool. Id prefer the duplicate_1 duplicate_2format but if it is much more complex then I am fine with the duplicate duplicate_1format.

@AudileF: I'm having a hard time coming up with a good reason why I would want the first duplicate to be changed as well.  Can you give us a bit more insight into the reason for this request?

Sure, this little Perl snippet should do it:
$ perl -pe 's/$/_$seen{$_}/ if ++$seen{$_}>1 and /^>/; ' file.fa
>1_uniqueGeneName
atgc
>1_anotherUniqueGeneName
atgc_2
>1_duplicateName
atgc_3
>1_duplicateName_2
atgc_4

Or, to make the changes in the original file, use -i:
perl -i.bak -pe 's/$/_$seen{$_}/ if ++$seen{$_}>1 and /^>/; ' file.fa

Note that the first occurrence of a duplicate name isn't changed, the second will become _2, the third _3 etc.
Explanation

perl -pe : print each input line after applying the script given by -e to it.
++$seen{$_}>1 : increment the current value stored in the hash %seen for this line ($_) by 1 and compare it to 1.
s/$/_$seen{$_}/ if ++$seen{$_}>1 and /^>/ : if the current line starts with a > and the  value stored in the hash %seen for this line is greater than 1 (if this isn't the first time we see this line), replace the end of the line ($) with a _ and the current value in the hash

Alternatively, here's the same idea in awk:
$ awk '(/^>/ && s[$0]++){$0=$0"_"s[$0]}1;' file.fa
>1_uniqueGeneName
atgc
>1_anotherUniqueGeneName
atgc
>1_duplicateName
atgc
>1_duplicateName_2
atgc

To make the changes in the original file (assuming you are using GNU awk which is the default on most Linux versions), use -i inplace:
awk -iinplace '(/^>/ && s[$0]++){$0=$0"_"s[$0]}1;' file.fa

Explanation
In awk, the special variable $0 is the current line.

(/^>/ && s[$0]++) : if this line starts with a > and incrementing the value stored in the array s for this line by 1 evaluates to true (is greater than 0).
$0=$0"_"s[$0] :  make the current line be itself with a _ and the value from s appended.
1; : this is just shorthand for "print this line". If an expression evaluates to true, awk will print the current line. Since 1 is always true, this will print every line.

If you want all of the duplicates to be marked, you need to read the file twice. Once to collect the names and a second to mark them:
$ awk '{
    if (NR==FNR){
        if(/^>/){
            s[$0]++
        }
        next;
    }
    if(/^>/){
        k[$0]++;
        if(s[$0]>1){
            $0=$0"_"k[$0]
        }
    }
    print
}' file.fa file.fa
>1_uniqueGeneName
atgc
>1_anotherUniqueGeneName
atgc
>1_duplicateName_1
atgc
>1_duplicateName_2
atgc

IMPORTANT: note that all of these approaches assume you don't already have sequence names ending with _N where N is a number. If your input file has 2 sequences called foo and one called foo_2, then you will end up with two foo_2:
$ cat test.fa
>foo_2
actg
>foo
actg
>foo
actg
$ perl -pe 's/$/_$seen{$_}/ if ++$seen{$_}>1 and /^>/; ' test.fa
>foo_2
actg
>foo
actg
>foo_2
actg

If this can be an issue for you, use one of the more sophisticated approaches suggested by the other answers.

In the awk example I worry that duplicate IDs could be created if you have contigs like >foo and >foo_1 in the original file. foo_1 wouldn't be in s, but might already have been printed out if foo had already been encountered.

@DevonRyan ah, yes that's a very valid point. Not just the awk one either, all approaches here break if a duplicate sequence already exists with a _N. In that case, the extra verbiage added by python is actually worth it :P. I added a disclaimer, thanks.

And here’s a solution using R (with the Bioconductor):
fa = ShortRead::readFasta(infile)
ids = as.character(ShortRead::id(fa))
fa@id = Biostrings::BStringSet(make.unique(ids, sep = '_'))
ShortRead::writeFasta(fa, outfile)

Sure, using biopython:
#!/usr/bin/env python
from Bio import SeqIO

records = set()
of = open("output.fa", "w")
for record in SeqIO.parse("foo.fa", "fasta"):
    ID = record.id
    num = 1
    while ID in records:
        ID = "{}_{}".format(record.id, num)
        num += 1
    records.add(ID)
    record.id = ID
    record.name = ID
    record.description = ID
    SeqIO.write(record, of, "fasta")
of.close()

Change output.fa and foo.fa. One doesn't need to explicitly change the .name and .description, but that's handy to prevent the original ID from not appearing still (after a space).
Explanation

records = set() This will create a lookup table of all IDs written
for record in SeqIO.parse("foo.fa", "fasta"): Open foo.fa as a fasta file and iterate over the entries in it. These are objects with an ID, name, description, and sequence attribute (the name and description are the same as the ID if not present).
ID = record.id Memoize the entry ID (e.g., 1_duplicateName)
while ID in records: As long as the ID has already been seen, keep looping.
ID = "{}_{}".format(record.id, num) Start adding increasing numbers after the ID, such as 1_duplicateName_1 and 1_duplicateName_2. This will continue until the ID has not been seen.
records.add(ID) Add the unseen ID to the set.
record.id = ID Update the ID, the .name and .description are the same. If you don't do that, then you get output like >1_duplicateName_1 1_duplicateName. That's not really a problem, but it's excessive.
SeqIO.write(record, of, "fasta") Write the record to the output file in fasta format.

A benefit of biopython is the easy ability to change formats, so instead of fasta one could have substituted Genbank or another format if needed.


text:
Calculating the charge of a peptide computationally
I was wondering how I can calculate the charge of a protein peptide (e.g. "RKTTLVPNTQTASPR") computationally in R or another tool.

Have you done any search already? What have you found?

A quick google search turns up protcalc, which is able to give a nice pH-dependent table of peptide charges (yours ranges from 3.1 at pH 4 to 1.5 at pH 10). It's on sourceforge, so hopefully the source code (and maybe a publication) is somewhere in there. Granted, this isn't R, but it's not clear from your post how necessary that really is.
BTW, in general "best and most accurate" turns out to be highly subjective (at least the "best" part of that).

Actually somewhere I read that a method to calculate the charge was taking in mind the pkAs of the free amino acids and not the residues on the chain. That's why I wrote "best". I think that the pkA change when the a.a is on the chain. Maybe I should have mentioned it. Also R environment will be preferred.

I would be surprised if pH still didn't affect things, though I agree that using a calculation based on free amino acid charges is going to yield wrong values.

Taking in mind the pH is standard but pH value is known and can be found. The thing is if there is an algorithm to calculate the pkA of the non free a.a

@J.Doe asking for "the best" is not a good question for this site since "the best" is almost always going to be a matter of opinion, and/or depend on the specifics of the sequence you are analyzing. See https://bioinformatics.stackexchange.com/help/dont-ask

Another quick Google search reveals the Peptides R package.
library(Peptides)
s <- "RKTTLVPNTQTASPR"
charge(s)

[1] 2.997683

Optional arguments to charge() are pH (default = 7) and pKscale (choice of nine, default = "Lehninger"). See ?charge for details.
Another (non-R) option is the EMBOSS suite program iep. I highly recommend installing EMBOSS, it provides a suite of useful command-line tools with output that can be easily read and processed by other sofware, including R.

I've used PropKa in the past to get isoelectric points. Pretty simple to use:
https://github.com/jensengroup/propka-3.1


text:
How to simulate "base error rate" in art_illumina?
I'd like to simulate 10% sequencing error using art_illumina. The simulator doesn't have a parameter that I can just give the 10%, but it has this:
-qs  --qShift   the amount to shift every first-read quality score by
-qs2 --qShift2  the amount to shift every second-read quality score by
                NOTE: For -qs/-qs2 option, a positive number will shift up
                quality scores (the max is 93) that reduce substitution sequencing
                errors and a negative number will shift down quality scores that
                increase sequencing errors. If shifting scores by x, the error
                rate will be 1/(10^(x/10)) of the default profile.

Heng Li's wgsim has an option for base error rate. Can I do the same for art_illumina?
Q: Does the "error rate" in the description mean sequencing error? If I set both qShift and qShift2 to 10, 1/(10^(x/10)) == 0.1, where x is 10, does that mean 10 is the value I should give to the simulator?

"a positive number will shift up quality scores that reduce substitution sequencing errors" -- I suspect you want to do the opposite: "a negative number will shift down quality scores that increase sequencing errors"

You need to specify an error profile to do this. However, you can shift quality scores of an existing profiles to get 10% error, which is easier than defining a complete profile.
For example, Illumina HiSeq 1000 has the error rate of 0.1% and you want 100 time more errors, so you should shift down the quality scores by 20 by setting both qShift and qShift2 to -20. The following code will simulate paired end reads with 10% error.
art_illumina -ss HS10 -i reference.fa -qs -20 -qs2 -20 -sam -p -l 100 -f 20 -m 200 -s 10 -o output


text:
Genome databases of microbial eukaryotes
I am interested in collecting the reference genomes of marine protists and would like to know which databases can be used for such tasks, and is there a single, centralized database of protistological genomic data?

Hi and welcome to the site. Please not that asking for "the most commonly used" isn't a good question for this site since it will almost always require opinion-based answers ("I use database foo!"). A better question would be to tell us what you want to do and ask how to do it instead of asking for the most popular DB. That said, please [edit] your question and clarify what you want to do. Are you looking for the reference genome of specific species? What kind of analysis are you thinking of? What species? Protists? Fungi? Algae? What kind of info would a database have?

You can find many protist genomes on ensembl
http://protists.ensembl.org/index.html


text:
How to isolate genes from whole genomes for phylogenetic tree analysis?
I have 446 whole Klebsiella Pneumoniae genomes I want to build a phylogenetic tree from. After reading about constructing phylogenetic trees it seems the only option for large numbers of genomes is to isolate a gene with low variability from generation to generation and use this gene to build a tree. For example Lars Jensen recommends using "16S rRNA [or] all ribosomal-protein-coding genes" https://www.biostars.org/p/1930/. What program isolates these genes of interest from the whole genome fasta files and can put them into a multiple alignment file? Or outputs them in a formate ready for a multiple alignment program such as Muave? The reason I say a multiple alignment file is because this it the type of file most phylogenetic tree programs take (I.E. clonalframe).

Hi Daniel, I would be cautious to build a phylogeny on 16s rRNA sequences as informative variability is essentially nil. I would suggest a good starting point being concatenated MLST genes, and work from there. You have the whole genome seems wrong to focus on the least variable gene.

Yes. Using 16S makes sense for distant species but not for closely related ones. it is extremely likely that there will be next to no variability at all.

If you haven't yet annotated the genomes, do so with Prokka as suggested in some of the other answers. Once you've got that, you can easily pull gene features out of a genbank by name using biopython, once you've identified a few sequences you want to use as the basis for your sequence typing.

There are a lot of ways to do this. I suggest using Prokka/Roary to produce a core genome alignment. There's a useful tutorial on the Roary website:
for file in *
do
    prokka --kingdom Bacteria --outdir "${file%%.*}"  --genus Listeria --locustag "${file%%.*}" "$file"
    mv "${file%%.*}"/PROKKA_07052017.gff GFF/"${file%%.*}".gff # use current date
done
roary -f Alignment -e -n -v GFF/*.gff

Alignment/core_gene_alignment.aln can be used as input for phylogenetic analyses

Roary is my go to but you have no say in what come out in the alignment. If you want to select your gene content use blast as I suggest. Concatenate the output and use an alignment software.

@terdon I added a code snippet to do the analysis. It can be used to isolate genes (plural). If the OP just wants to build a phylogenetic tree, there's no need to isolate specific genes. With almost 500 genomes from a single species there's unlikely to me much resolution from analyses of a single gene or from a handful of MLST genes.

Extract desired gene sequences using standalone blast
Simply provide a reference database with your desired output.
Set up your command and away you go. You can set the search up with a for loop for a batch of sequences. Command may look like
for f in *.fasta; do
   f=$(basename $f .fasta)
   blastn \
   -outfmt "6 sseqid qseq %" \
   -query $f.fasta \
   -subject reference.fna \
   > out/$f.fas
done

Watch the output as blast will spit out the gene as detected in + or - sense. If you want to only gather positive sense use -strandoption.
The default output I have here is tab output which requires a few sed commands to make into fasta.
sed -i \
-e 's/\s*$//g' \
-e 's/^/>/g'  \
-e 's/\s\+/\n/g' \
*.fas

Online alignment servers are an easy way to align small datasets e.g. EBI

It might be worth mentioning that this approach work for species with no splicing (so it should be fine for the OP), but isn't very useful for eukaryotic organisms. Also, why would you run multiple blasts in a loop like this? Why not give a multifasta input file to blast and run it once?

In my experience the blast out uses the gene names from the reference. and does not append the original genome name. I go about a few QC/ parsing bits and then combine the sequences.

I'm not sure what my "desired gene sequences" are supposed to be. What is my query?

Your query is the genome sequences you have $f will insert the name for each in a loop. your desired gene sequences will be in your reference file. This could be the 16s or MLST gene sequences. Whatever you want to build a phylogeny of.

My question is more how can I find those MLST gene sequences in the first place?

A quick google brought this up http://bigsdb.pasteur.fr/klebsiella/klebsiella.html . Have a look about there is usually a file full of genes and alleles. In my experience one (allele) set of the 7 genes should suffice.

@DanielHarris take one MLST gene and use that as a query to find the rest.

@AudileF blast can take a multifasta file as a query and the results for each sequence in the input file will be shown in a separate Query= section. The names are indeed the names taken from the input file, but that shouldn't be a problem as long as they're all unique.

So use one of these MLST genes such as from this link http://bigsdb.pasteur.fr/perl/bigsdb/bigsdb.pl?db=pubmlst_klebsiella_seqdef_public&page=downloadAlleles as a query. Two questions now... Do I use the entire locus file, such as gapA has many >gapA_#'s or just one >gapA_1? Second, how does the output of blast then be used to build a phylogenetic tree? Blast output is far from multiple alignment format I think.

@terdon must admit I never know. Does this work in the case of contig assemblies?

@DanielHarris 1) use one copy of each of the 7 genes. 2) you are correct the output can be quite jumbled. I usually use the -strandoption to gather the + and - sense sequences. Split them by taxa name. Order sequences in each file by name. remove fasta headers and space. Create new fasta header of file name. Then concatenate them all. For simplicity you could use a online msa server like on ebi.  Long winded and confusing. Im not a master computational biologist terdon might have a better suggestion.


text:
Which sequence alignment tools support codon alignment?
Sometimes it useful to perform a nucleotide protein coding gene sequence alignment based on codons, not on individual nucleotides. For example for further codon model analysis it is important to have full codons.
A widely used approach here is to perform a protein sequence alignment first and then impose this alignment to the nucleotide sequences using PAL2NAL, CodonAlign or something similar.
This is how transAlign or GUIDANCE (in codon mode) work.
The problem here is that you are discarding part of the information which could be potentially used for the sequence alignment. E.g. if you have slowly evolving low-complexity region adjacent to a quickly evolving one, the amino acid induced alignment could be wrong, while incorporating nucleotide sequence potentially allows to make the alignment more accurate.
I'm aware of two programs which can do true codon alignment. First, PRANK has a dedicated codon model, but it is rather slow and using it is overkill for certain problems. Second, Sequence Manipulation Suite can perform codon alignments, but only for a pair of sequences; also it's javascript based, therefore it is hard to run it for a large number of sequences.
Can you recommend any software for multiple codon sequence alignment? Preferably available for offline use.

Why not translate the sequence and perform amino acid alignment? True, you “discard” the codon isoacceptor identity but that can be mapped back trivially if needed.

As I write in the question, this way you loose some information. E.g. you can produce an incorrect amino acid alignment for low-complexity region, which you otherwise would be able to reconstruct MSA properly.
In principle I agree that this doesn't happen very often.

There's blastx and tblastx but I don't think they do true codon alignment.

@IakovDavydov how does a codon alignment 'loose some information' over a simple DNA-DNA alignment? It's easier to find homologous residues using proteins (instead of DNA), and using that information to align DNA triplets can only be beneficial.

@PejoPhylo I was not claiming that codon alignment looses some information. I was saying that when converting codons to amino acid some information is lost.

@IakovDavydov I'm sorry, but then I don't understand what you originally asked in this question. Couldn't you use PAL2NAL for your (multiple) codon alignment? It supports offline use

@PejoPhylo I don't think I can make it much clearer here than in the original question. But briefly: with PAL2NAL you loose some information, which might be important for some alignments. Most of the time that's not important, but if your dN is comparable with dS and GC-content stays constant, that could make a difference. This paper provides some details: https://doi.org/10.1093/molbev/msr272

@IakovDavydov thanks, I'll give it a read.

@IakovDavydov I don't think that your original question is correct.  If one is interested in the functional aspect (the protein), then any codon based aligner will work, for example MEGA.  By definition, the amino acids are being aligned. The question about the low complexity, slowly evolving region is mute. If you are interested in the codon differences and drift, or evolution of the codons, then yes - but, I think that using nt aligned by codon is done more from the perspective of bundling the nt into triplets.

I don't know of any transcript-to-transcript aligners that are able to do this, but LAST can align transcript queries to protein reference sequences using a specified frameshift cost. Here's the specific documentation for that option:

-F COST
Align DNA queries to protein reference sequences, using the specified
  frameshift cost. A value of 15 seems to be reasonable. (As a special
  case, -F0 means DNA-versus-protein alignment without frameshifts,
  which is faster.) The output looks like this:
a score=108 s prot 2  40 + 649
FLLQAVKLQDP-STPHQIVPSP-VSDLIATHTLCPRMKYQDD s dna  8 117 + 999
FFLQ-IKLWDP\STPH*IVSSP/PSDLISAHTLCPRMKSQDN

The \ indicates a forward shift by one nucleotide, and the / indicates
  a reverse shift by one nucleotide. The * indicates a stop codon. The
  same alignment in tabular format looks like this:
108 prot 2 40 + 649 dna 8 117 + 999 4,1:0,6,0:1,10,0:-1,19

The "-1" indicates the reverse frameshift.

I sent an email to the LAST mailing list about adding a frameshift penalty for transcript-to-transcript matching; I've been pleasantly surprised with the requested features that Martin Frith has added to LAST in the past. Unfortunately, in this case the problem is too difficult to sort out due to all the possible combinations that could happen, so it's unlikely to be implemented in LAST in the forseeable future (unless someone else writes that code).

Try MACSE v2 (https://academic.oup.com/mbe/article/35/10/2582/5079334) will align multiple protein-coding nucleotide sequences based on their amino acid translation while allowing for the occurrence of frameshifts

Thanks @user90!

I had no idea! thanks, the 'new' singularity provided here  works like charm.

Virulign does this for virus sequences, the publication is available here and github address here

Hi @kristoftheys welcome to the site, congratulations on Bioinfo paper. I've edited the post to include publication and github address.

@Michael many thanks for the completion. sorry for the inconvenience, next time better :)


text:
Human body sites knowledge base
I'm looking for a lightweight knowledgebase describing the human body to annotate disease sites.
I do not need a great level of detail, I just need a kind of basic organs/sub-organs taxonomy. I checked out the Foundational Model of Anatomy ontology, but it's over-kill for my purposes and really cumbersome (Protege fails when it tries to load FMA.owl).
Is there a simple online knowlegebase available for this?

I did manage to open FMA.owl in Protégé, by the way. It took a while and a considerable amount of CPU and RAM but it did work eventually.

Ok! Probably I have no enough HW resources to load FMA and then extracting a small subset from it. That's the reason why I was looking for something smaller.

I suggest you post a new question asking how to extract a subset of an existing ontology. Give us a link to FMA, explain what you need (for instance, how deep you want to go in the ontology's branches) and ask for help making a "slim" ontology from it.

I am the developer of Uberon and I would be happy to help you with what you need, either from Uberon, or from the FMA.
You mentioned you need something simpler than FMA. There are a variety of tools for creating custom subsets. Additionally, some ontologies provide ready-made subsets for particular purposes.
For example, we make a subset called 'basic.obo' from Uberon:
http://uberon.github.io/downloads.html#subsets
You can examine this directly in either oboedit or Protege5 (note Protege5 parses obo so you don't need to convert).
From your question, you didn't indicate you needed programmatic access, but if you do there are a variety of options, as well as more programmer-friendly JSON exports.
Most ontologies have a number of different options for web-browsing, you can see some of the options for Uberon here:
http://obofoundry.org/ontology/uberon.html
You mentioned annotating disease sites. There may be a tools (either desktop or web-based) for helping with this - e.g. providing autocomplete over a desired subset of an ontology.
You may also want to ask your question on a mailing list like obo-discuss. Unfortunately many ontology developers don't read stackexchange, but most are always keen to hear of community requirements, to help them provide better and more useful products.

Search in bioportal: For example:  SNOMED, UBERON.
Please, note that FMA is also available as RDF, csv , etc...
You could also try to convert FMA to obo and try to open it with OBO-edit

SNOMEDCT seems to be a full ontology for any type of clinical data, not only anatomical (e.g. it includes terms such as "smoker" or "observation") which would make it far more complex than the FMA the OP tried. Likewise, Oberon seems to be a cross-species ontology which, again would be more, not less complex than FMA. If I understand correctly, the OP wants something less complete and less complex than FMA.

Yes @terdon, you are right. I'm looking for something simpler and smaller than FMA.

@Pierre, I'm afraid that Oboedit performs even worse than Protege


text:
using python to write bioinformatics pipelines tutorial
I was wondering if there is a tutorial or a small code snippet to understand how to write bioinformatics pipeline using python, for example

use a aligner (say hisat)
get the output and process it using samtools

I was able to use subprocess from python2.7 for this purpose using samtools but i am not able to link both the processes.i.e given path(which I can use argparse) for directory with fastq files output would be processed bam.
sample code for samtools sam to bam :
import subprocess
subprocess.run(['samtools','view', '-bS',
                '../some.sam',
                '>',
                '../some.bam'])

If you ever need more elaborate pipelines in python, then I would strongly encourage you to look into snakemake, which is pretty commonly used for this purpose.

I agree that snakemake is worth learning. If you want to use subprocess, I found the following tutorial just yesterday, which looked quite clear: http://crashcourse.housegordon.org/python-subprocess.html

@bli this was useful

@DevonRyan I looked snakemake it is available only in python3 is that correct?

@novicebioinforesearcher Yes, that's correct.

@novicebioinforesearcher Note that if you need to program in python2, you can still do it in external scripts and call them from snakemake like you would do for any other program.

@bli how can I do that, sorry I am still learning stuff here could you share an example

@novicebioinforesearcher These comments are not adapted to give examples. What I mean is that you will still need python3 for snakemake, but if there are things you can only do in python2, you can write whatever computation you want to do in python2 in the form of an independent program (presumably accepting command-line arguments). Then, you will call this program from within snakemake (or whatever workflow management tool you happen to use, or a python script using subprocess) exactly like you would do for samtools, hisat or any other program.

I'm curious as to why you are averse to writing your pipeline in bash? @novicebioinforesearcher

@quantik because I am learning to write codes in python and also because (i may be wrong here) in case you have a greater number of samples bash would not be a good way to work(foreg 40 samples to be processed) i can use os.walk etc in python to submit jobs correct me if i am wrong

I agree that using a specialized tool is probably a good idea.
Nevertheless if you want to stick with Python, I suggest using plumbum instead of subprocess. It has a very nice syntax for this kind of problems.
from plumbum import local

# load a command
samtools = local['samtools']

# you can then easily redirect
(samtools['view', '-bS', '../some.sam'] > '../some.bam')()

# it is also easy to pipe
hisat2 = local['hisat2']
chain = hisat2['-x', 'genome.index', '-1', 'reads_R1.fq.gz', '-2', 'reads_R2.fq.gz'] | samtools['sort'] > 'reads_vs_genome.bam'

chain()

Taking a different tack from other answers, there's lots of tools for pipelines in Python. Note: there was a time when people would use "pipeline" to refer to a shell script. I'm talking about something more sophisticated that helps you decompose an analysis into parts and  runs it robustly.

Snakemake is my favourite. It's (nearly) pure Python and can generate reports.
Nextflow is growing in popularity and is pretty straightforward
Ruffus used to be reasonably popular, seemed fine when I used it
Bpipe is for bioinformatics
Airflow is a more industrial "big" solution but in wide use
See this big list of pipeline systems

I'm sure you can find something there.

BioPython has some good tools for processing reads and alignments.
http://biopython.org/DIST/docs/tutorial/Tutorial.html
There is a python library wrapping samtools so many of the samtools calls can be used directly as python objects and calls
https://pysam.readthedocs.io/en/latest/
I would use subprocess to call the aligner and specify the output to a bam file that you have named and then read that bam file with pysam to do the analysis that you are interested in.
If you provide more specifics about the analysis or aligners I can help you more but you are in the right track just need to connect the two calls together.
example:
import subprocess
import pysam
import os

for root, dirs, files in os.walk(rootPath):
    for filename in fnmatch.filter(files, ".fastq"):

        subprocess.run([`hisat`, `-f`, filename `-o`, `some.bam`])

        #get depth
        depth = pysam.depth('some.bam')

edited code snippet based on OP comments about executing on multiple fastq files

If you're just doing alignment and conversion to a sorted BAM file, there's no need to run it through python. A simple pipe on the unix command line works just as well (and probably runs faster):
hisat2 -x genome.index -1 reads_R1.fq.gz -2 reads_R2.fq.gz |
  samtools sort > reads_vs_genome.bam

thank you @gringer, I am familiar with unix pipe. I am trying to learn to create a pipeline where when I give path of a directory which has fastqs I intend to submit .sh files to cluster and get the desired results . for example python my_pipeline.py -dir /pathtofastq  would submit jobs to cluster here i will use bunch of argparse and subprocess and predefinec bsub commands.

Look into python's os.path.walk.  This will allow you to take in the dir and then iterate over all the fastq files in the directory.  http://www.pythonforbeginners.com/fnmatch/os-walk-and-fnmatch-in-python

Have a look at luigi for writing pipelines in python. You can found BioLuigi if you like :-P

could you please share link for BioLuigi I could not find it online

You need to start it ;-)

haha good one!!


text:
Convert Reactome Protein IDs to Pathway IDs?
I have a list of about 50k 'Protein IDs' from Reactome. Is there a simple way to get all the corresponding 'Pathway IDs' for each protein? What is the best service to use? (I'm guessing I can use the Reactome API, but I don't necessarily want to hit that 50k times...).
Protein IDs from reactome look like this:
R-HSA-49155
R-HSA-199420

The corresponding Reactome 'Pathway IDs' for those Protein IDs would be:
R-HSA-49155  R-HSA-110331
R-HSA-49155  R-HSA-110330
R-HSA-49155  R-HSA-110357
R-HSA-199420 R-HSA-1660499
R-HSA-199420 R-HSA-202424
R-HSA-199420 R-HSA-199418

Please [edit] your question and give us a few IDs as examples we can use to test our approaches. And why don't you want to use the API? This sort of thing is precisely what APIs are for. You can always limit it to only N requests per minute or whatever they need.

Yes there is a way in R to do this with the "reactome.db" library, but it is unclear from your question what these protein IDs are, and what organism you use. You'll probably have to convert you protein IDs to eg first, if you use human data you'll need for example "org.Hs.eg.db" library.

I must say that I tried to answer this question using reactome.db of Bioconductor and I couldn't. The current version (1.59.1) doesn't provide the reactome ID for a given ENTREZ ID, only the pathways IDs. I submitted a question to the support site

This should be an easy task for biomaRt in R, hoewever, I found out that they use the reactome_gene IDs as if they were pathways IDs. So it doesn't work properly.

Thanks for trying... I guess the API is the way to go? I just don't like issuing 50k web-requests, although I guess it isn't really a problem. Just created a bounty for the lulz.

If you don't mind hitting it 50k times and are OK with python3...
from urllib import request
import json

def getPathways(proteinID):
    baseURL = 'http://reactome.org/ContentService/data/query'
    PathwayIDs = set()
    try:
        response = request.urlopen('{}/{}'.format(baseURL, proteinID)).read().decode()
        data = json.loads(response)
        if 'consumedByEvent' in data:
            for event in data['consumedByEvent']:
                PathwayIDs.add(event['stId'])
        if 'producedByEvent' in data:
            for event in data['producedByEvent']:
                PathwayIDs.add(event['stId'])
    except:
        pass
    return PathwayIDs

Usage would then be something like:
l = ['R-HSA-49155', 'R-HSA-199420', '']
for rid in l:
    ids = getPathways(rid)
    for _ in ids:
        print("{}\t{}".format(rid, _))

Which would produce:
R-HSA-49155 R-HSA-110239
R-HSA-49155 R-HSA-110240
R-HSA-49155 R-HSA-110238
R-HSA-49155 R-HSA-110356
R-HSA-199420    R-HSA-8948800
R-HSA-199420    R-HSA-6807106
R-HSA-199420    R-HSA-6807206
R-HSA-199420    R-HSA-6807126
R-HSA-199420    R-HSA-8847968
R-HSA-199420    R-HSA-8850997
R-HSA-199420    R-HSA-2321904
R-HSA-199420    R-HSA-8948775
R-HSA-199420    R-HSA-8944497
R-HSA-199420    R-HSA-6807134
R-HSA-199420    R-HSA-8850945
R-HSA-199420    R-HSA-8873946

Note that this will silently ignore invalid or missing IDs such as '' (that's the try and except above. Note also that these are different pathway IDs than what you provided in your example. The main reason is that the protein IDs you showed are not always involved in the pathways IDs you showed (in my example, they always are).

This looks great, but for some reason fails here: http://plantreactome.gramene.org/ContentService/ when I query with this ID: R-OPU-1129769. However, if I use that ID I do get lists of pathways here: http://plantreactome.gramene.org/ContentService/data/pathways/low/entity/R-OPU-1129769?speciesId=9007655
Sorry for using human reactome as an example when I wanted data from plant reactome (I assumed it would work equally well in both places).

On Reactome website they have at the download page mapping files (uniprot, ensembl, etc.), but unfortunately not for the protein IDs you are using (stable identifiers).
I had contact with their helpdesk, and they sent me a file containing all protein IDs to the pathways. Exactly what you need. I have asked them if they wanted to put it on their download page as well, but not sure if they want to do this. Meanwhile you can ask for the file as well, or get it from my google drive.
I assume you know how to get your Protein IDs and pathways from this file, using e.g., R?


text:
scanBam from Rsamtools is not importing one of my reads into R
I have this read in my BAM file. It maps on chromosome 1.
I open this BAM file in IGV, and I can see the alignment on chromosome 1.
But when I open this file in R with Rsamtools:
bamContigsCel <- Rsamtools::scanBam('output/alignment/pacbio/bwa/ref     /bristolAssemblySorted.bam', param = Rsamtools::ScanBamParam(what = Rsamtools::scanBamWhat(), flag = Rsamtools::scanBamFlag(isMinusStrand = FALSE), tag = bamTags))[[1]]

I then check if the read maps to chromosome 1 in my R object but I cannot find it.
bamContigsCel$rname[bamContigsCel$qname == '000000F|arrow']
[1] II II IV
Levels: I II III IV MtDNA V X

But if I look at the BAM file, it's there. Why isn't Rsamtools importing my read into R?
$ samtools view bristolAssemblySorted.bam | grep -n '000000F' | head  -c80
000000F|arrow   2064    I   336331  7   2926310H260M1774267H    *   0 0   GAAGCTGTCTAAACTTTGGC

Cross-posted on biostars

Careful, your code is using an almost certainly wrong path to the BAM file.

Note the flag; that read is mapped in a reverse-complemented manner, so isMinusStrand = FALSE is filtering it out. I tested this by making a BAM file with only that read:
@SQ SN:I    LN:1000000
000000F|arrow   2064    I   336331  7   2926310H260M1774267H    *   00  *   *

Then in R:
> library(Rsamtools)
> blah = scanBam("foo.bam", param=ScanBamParam(what=scanBamWhat(), flag=scanBamFlag(isMinusStrand = TRUE)))[[1]]
> blah$qname
[1] "000000F|arrow"
> blah = scanBam("foo.bam", param=ScanBamParam(what=scanBamWhat(), flag=scanBamFlag(isMinusStrand = FALSE)))[[1]]
> blah$qname
character(0)

Thanks! I don't know why I was filtering out reverse-complemented alignments! Can I just use bam <- Rsamtools::scanBam('output/alignment/pacbio/bwa/ref/bristolAssemblySorted.bam')[[1]] to get all the information in the BAM file?

Yup, that should work!


text:
COG Annotation - Dealing with genes assigned to two or more COG categories
I am dealing with a list of genes which have been selected from a gene enrichment analysis. In order to see what kind of genes are overrepresented, I ran eggnog-mapper to do an orthology assignment of each gene against a bacterial database (provided in eggnog).
After running eggnog you get a .csv file containing for each gene which is the closest ortholog, and a list of GO, KEGG and COG terms associated.
From there, I would like to do a Fisher-test (also considering Bonferroni correction) to highlight COG categories overrepresented in our dataset compared to our reference genome.
Finally my question: How do you deal with genes assigned to multiple COG categories (same would apply for GO terms). You count +1 to each of the COG categories assigned? Do you just choose one category per gene?

You add 1 to the count of each COG category. You are looking for over represented COG categories so you must count them all. Many genes will be assigned to multiple categories. In fact, a majority of human genes are given multiple GO cellular function annotations. I have worked considerably in this field and one of the surprises I found is how common protein multifunctionality is. But I digress.
The point here is that you want to see whether any COG categories are over represented in your test set. Whether a gene is annotated to one or to a hundred categories is irrelevant. You are counting genes per COG, so you need to count each of the COGs assigned to a gene separately.


text:
Y Chromosome Aligned Reads in scATAC-seq data from a female-derived cell line?
I'm working with scATAC-Seq data on the K562 cell line, which is supposed to be derived from a female patient. While following the scATAC-seq data analysis pipeline, after performing bowtie alignment they recommend filtering out all reads aligned to mitochondrial DNA and on the Y chromosome.
Out of curiosity, I decided to count how many reads aligned to chrY and found that it can be quite high - chrY has approximately 10% as many reads as a similarly-sized chromosome (chr19 in hg19).
Is it normal to find reads mapping to chrY when the cell line is female-derived, or did I mess up somehow?

There are homologous regions between X an Y chromosomes: https://en.wikipedia.org/wiki/Pseudoautosomal_region
It is therefore normal to have some female-derived reads mapping in Y chromosome.
You should probably check what proportion of such reads fall in other parts of the Y chromosome than pseudoautosomal regions.

This is it -- I just noticed that most prominent peak corresponds to PRKY, a peudogene in the pseudoautosomal region that's similar to PRKX. Not all of my reads in chrY were pseudoautosomal, but all of my MACS2-called peaks were. Thanks!


text:
Software to produce a table of post-translational modifications from a peptide list
Does anyone know if there is a program/library/script in R or Python that takes as input a list of proteins/peptides and a list of post-translational modifications (PTMs; like oxidation of methionine and acetylation of cysteine), and as output returns a matrix/list with the imported peptides and all possible locations of the selected modifications?

Why would you need that? Maybe you can work with Uniprot modifications listed at their webpage instead of creating "new" ones?

I don't know of a tool that can do exactly what you want. However, there are various tools that can predict specific classes of post translational modifications. For example:

Phosphorylation: GPS

Computational prediction of phosphorylation sites with their cognate protein kinases (PKs) is greatly helpful for further experimental design. Although ~10 online predictors were developed, the PK classification and control of false positive rate (FPR) were not well addressed. Here we adopted a well-established rule to classify PKs into a hierarchical structure with four levels.

Mitochondrial import: Mitoprot

MitoProt calculates the N-terminal protein region that can support a Mitochondrial Targeting Sequence and the cleavage site. A complete description of the method to make the prediction is available in: M.G. Claros, P. Vincens. Computational method to predict mitochondrially imported proteins and their targeting sequences. Eur. J. Biochem. 241, 779-786 (1996).

Acetylation: "PAIL: Prediction of Acetylation on Internal Lysines"

In this work, we present a novel online predictor for protein acetylation sites prediction of PAIL, Prediction of Acetylation on Internal Lysines. We have manually mined scientific literature to collect 249 experimentally verified acetylation sites of 92 distinct proteins. Then the BDM (Bayesian Discriminant Method) algorithm has been employed.

You can use tools like the above and then combine the results to get the sort of matrix you are asking for.


text:
How to interpret contig-alignment.psa produced by velvet
I'm using velvet to align given reads of RNA to given CDSs (i.e. coding areas and genes) of an organism, so I can generate gene-expression profiles. But after using velvetg out-dir/ -alignment yes, velvet produces contig-alignment.psa which is a file in a strange format.
The file for every contig contains zero or more records of 7 integers, and I cannot understand their meaning.
The following is an excerpt of the file, and you can see the entire file here:
>contig_1
1    193    1939    1939    6   194 2
195 244 1   1   6   22  71
>contig_2
1   84  1   1   6   59  142
86  170 1935    1935    6   22  106
172 285 1935    1935    6   108 221
>contig_3
1   98  2   2   6   1   98
100 321 2   2   6   100 321
334 415 1204    1204    6   1   82
>contig_15
1   23  3   3   6   84  106

Is that the only file produced? The .psa format is not like what you show (I think, anyway, I haven't used it before but that's what I could find). Are there really no other output files? If not, it sounds as though that might be a bug.

Of course there are other files that also appear while using velvet without the option -alignment yes. But related to contig alignment I can only see this file

I don't think velvet is doing what you think it is. Velvet is a genome assembler. To quote the velvetg command help: "-alignments <yes|no> : export a summary of contig alignment to the reference sequences (default: no)" If you want to align RNA-seq reads to a transcriptome and get "gene-expression profiles" use a tool such as bowtie2 or Salmon.

Thank you, I used Last aligner, but the reason I wanted to use velvet was because It generates coverage ratio that could indirectly imply gene expression.

@hhoomn It does in a very roundabout way and it won't give any measure of relative expression between genes. I think you need to explain what exactly it is you want to do in terms of a biological question. It seems to me you're asking the wrong question, here.

The question is that what does this file contain, and the answer is yet unclear, but I solved what I needed to do with another tool, if you want to know what I was doing see here

It looks like you are trying to do a transcript isoform expression analysis. Velvet is a genome assembler, and really not designed for this task. The tools you are looking for are kallisto and salmon.


text:
How can I compute gene expression for a set of RNA reads?
I'm trying to compute a gene expression profile for an organism. I have gene nucleotide sequences of the mentioned organism stored in a fasta file and a set of paired reads stored in two separate files with the same fasta format. Now I want to compute the coverage of every gene by the given reads. That is, on average how many reads cover a given segment of the sequence.
Is there any tool for this specific purpose?
I'm trying to use the velvet sequencer but I have encountered a few problems.

It's unclear if your paired-end reads are actually in fasta format, I'll presume that they're in fastq instead.
The easiest tool to use is salmon, which nicely deals with things like multimapping. If you're trying to judge the quality of an assembly, then I recommend having a look at transrate, which uses some related methodology for assessing contig quality.

Yes they are in fasta format, but they can be easily transformed to fastq with  uniform quality.

What happened to the FASTQs? Technically you could convert to FASTQs, but if the reads are not FASTQs to begin with, it's often a red flag indicating other problems.

Thank you I used Last aligner for the work, and calculated coverage ratio manually (and considered that as proportional to gene expression).

In addition to salmon that was already mentioned, you can also try kallisto or RSEM, which are also fairly popular/respectable and will work with a transcriptome FASTA.

It's not possible to compute absolute expression from RNASeq reads if they are processed in the usual way, where a sequencer produces the same number of reads regardless of the input RNA amount. At best, RNASeq will give you an indication of proportional expression within a single sample. For this reason, relative expression (i.e. that used by differential expression tests) is easier to determine than absolute expression.
The closest approximation to absolute expression is to generate an expression relative to the average expression of a set of housekeeping genes, but I don't think there's a universal set that has been decided on. Gene expression, even for common housekeeping genes, can vary depending on the environmental conditions of the cell. For example, GAPDK is involved in immune cell activation.
However, as long as experimental conditions are similar, and you're not planning on looking for statistical significance, the proportional expression can still give qualitative insights into how cell populations behave in relation to other populations. DESeq2 provides a variance-stabilising transformation function that minimises variation for small-count genes, assuming that each sample has roughly the same total expression. I have found that I get better outcomes / comparisons from this transformation when carrying out a further adjustment to account for gene length (i.e. divide by the length of the longest transcript for each gene). See our Th2 paper, section "Read mapping and differential expression analysis" for more information. The "transcripts-per-million" values produced by Kallisto and Salmon provide measures similar to this.
If, on the other hand, you were able to modify the experimental design, single cell sequencing (or "known cell count" sequencing) can be used for determining absolute expression: use a spiked-in transcript that is added in proportion to the cell count, so that results can be compared in proportion to the expression for that transcript.


text:
Where can I find a database that has phenotype information together with associated SNPs?
I have been searching for a long time and the furthest I got is some database with the functional description for genes. Then, I have to parse these descriptions manually to figure out the association. Moreover, for some species, there are not even any functional descriptions.
I need this information for mice, rice, and Arabidopsis Thaliana. Can anyone help me?

What do you start with? Do you have a list of phenotypes you want to match to SNPs or a list of SNPs to match to phenotypes? Or do you have neither and need to find both?

@terdon I have both SNPs and phenotypes. I want to know the association between them. However, I don't want to do this with GWAS, because I need to go into wet-lab to verify these findings anyway. I wonder if there is a database with wet-lab verified associations already.

International Mouse Phenotyping Consortium is building a database of phenotypes and knock-outs of mouse. I believe that this database will be fairly complete (20000 knock-outs), but these are knock-outs, not SNPs... There are several mouse GWAS studies, but I am not aware of a database that would pull all the results together.
Arabidopsis big GWAS project is called 1001 genomes project. A lot more information about Arabidopsis from individual experiments is collected and maintained by The Arabidopsis Information Resource, however this resource is subscription based, therefore what you need might not be available for free.
As mentioned, SNPs of rice were identified and linked to different rice groups in The 3,000 rice genomes project. There is also GWAS done on 512 individuals, not sure if it is sufficient for reliable conclusions (even they claim that they found supersigificant links).
Is this what you were searching for?

A few days ago, I was trying to find some GWAS datasets to download. I hope this site for 3000 rice will help: http://snp-seek.irri.org/. You can click Genotype in the index page for "Query for SNPs from the 3000 genome project".
UPDATE: With the help of user manual, I can get IRIS ID from Genotype. According to this paper, iris id may associate with some phenotype information. But the website of IRIS is under maintenance, I can not get more information about phenotype.

but it is not linked to phenotype, no? When I checked they just associated it is different populations...

Hi l0o0, thanks for your answer and welcome to Bioinformatics Stack Exchange. I have edited your answer to improve the grammar slightly; if you can provide more information about snp-seek (e.g. whether or not it helps for finding associations with phenotypes), please edit your answer and add that in.

Have you looked at Ensembl? At least for mice, they have pretty rich data for SNPs e.g.
http://www.ensembl.org/Mus_musculus/Variation/Explore?db=core;r=2:3225281-3226281;v=rs27096498;vdb=variation;vf=802773
Assumes you have known dbSNP ids (or 'rs' numbers) for your SNPs. If not you could use their VEP tool to annotate them:
http://www.ensembl.org/Homo_sapiens/Tools/VEP?db=core
VEP is also available for rice and Arabidopsis:
http://plants.ensembl.org/Oryza_sativa/Tools/VEP?db=core

We also have variation and phenotype data for arabidopsis in Ensembl. The variation data is from the 1001 project and we took the phenotype data from Atwell et al., 2007 (http://europepmc.org/abstract/MED/20336072)


text:
Using data mining of papers in order to derive genomic connections
Are there any examples of researchers data mining articles and papers (e.g., from pubmed or google scholar) in order to derive relations between genes and diseases?

Do you have a specific question about text mining?

@jaslibra, it's a general question about this approach. For example to decide if a gene is tumor suppressor one. What kind of mutation are associated with pan-cancer, predict if a gene is oncogene. Basically, it has endless potential

This is really way too broad to be a good fit for this site. Unless you'd be satisfied with a yes, this can only be answered with a list fo examples. The simple truth is that everyone reads articles to learn the state of the art and then might move on to the lab to validate whatever inference they made. Please try and [edit] to make it more specific.

@terdon the question is whether a person reads it or a machine. That's the whole point.

The short answer is yes.  You could scrape pubmed for all the papers that are based on scraping data.  ;)

@0x90 both. Depends on what you need to do, that's why it's too broad. And you are not actually asking that, by the way. You are asking: Are there any examples of researchers data mining articles and papers [. . . ] in order to derive relations between genes and diseases? If you want to know how data mining works, that would be a very different question. If you want to know whether it is a machine or a human, you would have to give a specific example since both do it in different cases.

I'll provide an answer here, since the question got locked for whatever StackOverflow-ish reason. You might be interested in SmartFigures: http://smartfigures.net/ which uses data in figures to annotate an ontological database of genes, diseases, etc. See: https://vimeo.com/158555231 and https://vimeo.com/170011946 for more detail. Hope this helps and welcome to Stack Overflow!

If you feel this question was put on hold unfairly, please open a meta discussion. I am open to being corrected. I still feel this is very broad if you don't edit to limit to specific examples. I mean, of course researchers get disease associations from papers! Where else would they come from?

Also, the best way for you to reopen this question is to edit it and make it ask something specific instead of a vague "are there any examples". Explain exactly what you are thinking of (you made some more specific points in your comments, for example) do some research into what text mining is and how it can be applied and ask something specific. Or, you could change your question from "Are there any examples" to "how can this be done" which would also make it answerable.

There are many databases that have used publication scraping for oncogenic gene fusions. There are publications for the individual methods they used for their scraping and aggregation.

COSMIC - http://cancer.sanger.ac.uk/cosmic
TICdb - http://www.unav.es/genetica/TICdb/
ChimerDB - http://ercsb.ewha.ac.kr/fusiongene/ -
https://www.ncbi.nlm.nih.gov/pubmed/27899563

This existed as a closed silo, at least in 2015. Qiagen has a team of hired students and Post-Docs for collating research papers into their Knowledge Base, an extensive database that is integrated into a few of their commercial products. Qiagen's claim is that by providing a consistently-structured and well-formatted database, the process of research discovery is accelerated. I've got a few notes from a seminar they gave us a couple of years ago:

The curation is better than a free / public database because free databases typically have wrong or inconsistent data
Their knowledge base is an army of hired MDs and PhDs, 5.3M "findings" as of Nov 05 2015
Scientific articles have unstructured data, so expert knowledge is required to process the data into a "finding" for their knowledge base

I was already doubtful about their approach, and this doubt was magnified when they started getting into the details of their process. I started to lose trust and write my own thoughts about it:

The knowledge base they have created is private
The approach is to use a function (i.e. MDs/PhDs) to convert research articles into scientifically-validated associations. However, the function is effectively a black box; Qiagen allow the experts to decide for themselves how a particular article should be converted into a finding.
The findings in the knowledge base have no strength or likelihood.
The Experts were encouraged to only enter findings that were "generally accepted by the scientific community. but it seemed like only a single expert was required to get a finding into the knowledge base (i.e. findings weren't curated by other people).
When using their software (e.g. IPA), p-values of 10^-15 and 10^-22 are reported; these numbers are meaningless.
The team didn't look at methylation, phosphorylation, etc.; plans are in place to hire additional experts to look at this at a later date.

After listening to their seminar, I realised that Qiagen's idea of a general database that attempts to capture everything isn't going to work well. Experts will disagree about whether a particular association is valid, relevance (or correctness) will differ depending on the area of interest, and the established correctness of a piece of information can change over time.
It would be possible to create something that might work by adding a lot more metadata to all the "findings" in such a database, but it would need to be a huge, public, collaborative effort. Something on the scale of SRA/ENA, with encouragement from the scientific community to enter details into the database prior to publications being accepted/indexed. If the people who wrote papers entered their findings in such a database, other researchers wouldn't need to be hired to try to understand papers and guess at the findings of published research.

Thomson-Reuters is doing a similar thing, but I agree that "an association" is a complex thing and perhaps not to be taken at face value. Protein A may interact with protein B under lab conditions but whether this is a generally correct and useful thing is unclear.


text:
alternatives to MEDIPS to analyse MeDIP datasets
MEDIPS is an established tool with functions for the quality control and analysis of data derived from immunoprecipitation (IP)-seq samples, like Methylation IP sequencing datasets.
I would like to know if there are any other tools I should consider nowadays that could do a better job and hence replace MEDIPS for the analysis of MeDIP datasets.

I guess QSEA should be an obvious answer:

QSEA (quantitative sequencing enrichment analysis) was developed as
  the successor of the MEDIPS package for analyzing data derived from
  methylated DNA immunoprecipitation (MeDIP) experiments followed by
  sequencing (MeDIP-seq).

See http://bioconductor.org/packages/release/bioc/vignettes/qsea/inst/doc/qsea_tutorial.html


text:
Is Wikidata a reliable KB for genomic and drugs annotations?
I'm looking for a Linked Open Data approach to annotate a dataset with human genome and drugs information.
According to the Linked Open Data cloud, there are a lot of interconnected RDF vocabularies devoted to life sciences but, unfortunately, a lot of them seem out of date.
I recently found Wikidata:

Wikidata is a free and open knowledge base that can be read and edited
  by both humans and machines.

It seems exactly what Linked Data should be and it seems continuously updated. At a first sight, it covers very well genomic and drugs knowledge.
Does Anyone know if Wikidata is a reliable source for genomic and drugs annotation?

The question is what you mean with "reliable". For bacterial genome annotations there is one Wikidata based platform called WikiGenomes (publications).

Hi konrad, thanks for your answer and welcome to Bioinformatics Stack Exchange. If you would like your answer to collect more votes, it is helpful to add more information about any links that you provide. What does WikiGenomes provide in addition to annotations? How reliable do you consider it to be?

Thank you konrad. With "reliable" I mean regularly updated and curated from trusted sources. PS: I'm looking for human genome annotations. I'm going to specify that also in the question