text:
Duplicate long hits from PSI-BLAST
I had a protein Refseq ID and I PSI-BLASTed this sequence against Refseq database. We all know that the Refseq is a Reference sequence database and it shouldn't have redundancy. After BLASTing my sequence, at first iteration I got 1000 hits and among them there were a lot of redundant sequences!
My sequence had 241 amino acids and I found a lot of sequences with 100% identity, 100% cover and 0 E-value exactly the same as my sequence but with different IDs. All of these IDs were from RefSeq! In other iterations and after adjusting format options, I got this redundancy with other sequences from other species. My sequence is related to a chain of a multichain protein (E.coli fumarate reductase)
For example, when using NP_418578 as a query, I found WP_078165098.1, WP_064226696.1, WP_062863447.1, WP_001401474.1 and other that were identical.
I want to know what is wrong with Refseq. Is it really a Reference sequence database? If it is, where does this redundancy come from and why?

I assume you were using the protein version of the database, right? And were all the hits from the same species? Please [edit] your question and show us a specific example. Show us your query sequence and a few of the redundant results.

Yes, I am using the protein version of the database and my Refseq ID is: NP_418578.1  other Refseq IDs the same as this one are: WP_078165098.1, WP_064226696.1, WP_062863447.1, WP_001401474.1 and others

@Sara Maybe you could post a screenshot of your results (and queries). But having a short sequence of 241aa in several sequence might indicate good conservation (maybe it is a domain?)

@Sara Please [edit] your question to add extra information. Comments are easy to miss, hard to read and can be deleted without warning.

This is what is known as a feature, not a bug. Note that your identical proteins all have accessions starting with WP_. These are special "non-redundant proteins". Many sequences—particularly bacterial sequences—are identical between various different species so having a separate RefSeq entry for each of them would be inefficient. Therefore, RefSeq combines multiple such proteins into a single WP_ record. This is documented here (emphasis mine):

A new type of RefSeq protein record which represents non-redundant protein sequences was introduced in mid-2013. This record type was introduced to address a growing issue with redundancy in the Prokaryotic RefSeq protein dataset that coincided with a significant increase in bacterial genome submissions from individual isolates and closely related bacterial strains. For example, a large number of high-quality bacterial genomes may be submitted during a disease outbreak. The submitted sequences may reflect pathogen evolution during the course of the outbreak but the majority of the encoded proteins from these genomes may be identical to each other. As RefSeq includes these genomes, per community requests, this resulted in increased redundancy. By representing identical proteins using a single non-redundant protein accession number (with the prefix 'WP_'), redundancy in the database is significantly reduced.
[ . . . ]
Because a non-redundant protein sequence may be found in RefSeq genomes from multiple species, the organism information provided on the protein record reflects the lowest-common taxonomic node ranging from the genus species level to super-kingdom. A non-redundant protein record that provides organism information at the level of a genus, family, or even super-kingdom does not mean that the protein is found in all RefSeq genomes below that taxonomic classification. It only indicates that the protein is found in more than one genome of different species for which the genus, family, or super-kingdom classification is the lowest common taxonomic node.

So, your query was NP_418578.1, anaerobic fumarate reductase catalytic and NAD/flavoprotein subunit from E. coli strain K-12, substrain MG1655. The first thing to notice is how specific that is. This is the protein found from one specific substrain of one specific strain of one specific bacterial species. It is reasonable to expect that there will be identical sequences from many, many closely related species. Both from, most probably, all other strains and substrains of E. coli and from other, related bacteria.
Now, the specific sequences you mention are actually slightly different and not 100% identical. Below is a multiple alignment of NP_418578.1 and the 4 WP_ sequences you mentioned. Note that each of the 5 entries is slightly different. Each has one residue that differs from the rest. Look for the : in the identity line, there are 4 : and all others are * (I am only showing the relevant alignment blocks here and have removed those where all 4 sequences were identical):
WP_001401474.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
WP_062863447.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
WP_064226696.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
NP_418578.1         MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
WP_078165098.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAIAQDH
                    *******************************************************:****

[ . . . ]

WP_001401474.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
WP_062863447.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
WP_064226696.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHIRGLVAMNMMEGT
NP_418578.1         KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
WP_078165098.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
                    ***********************************************:************

[ . . . ]

WP_001401474.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
WP_062863447.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKISQAFWHEWR
WP_064226696.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
NP_418578.1         GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
WP_078165098.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
                    **************************************************:*********

WP_001401474.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGIDPVKEPIPVRPTAHYTMG
WP_062863447.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
WP_064226696.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
NP_418578.1         KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
WP_078165098.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
                    *****************************************:******************

[ . . . ]

Your sequence (NP_418578.1) is only identical to one WP_* multi-species sequence, WP_001192973:
WP_001192973.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
NP_418578.1         MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
                    ************************************************************

WP_001192973.1      DSFEYHFHDTVAGGDWLCEQDVVDYFVHHCPTEMTQLELWGCPWSRRPDGSVNVRRFGGM
NP_418578.1         DSFEYHFHDTVAGGDWLCEQDVVDYFVHHCPTEMTQLELWGCPWSRRPDGSVNVRRFGGM
                    ************************************************************

WP_001192973.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
NP_418578.1         KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
                    ************************************************************

WP_001192973.1      LVQIRANAVVMATGGAGRVYRYNTNGGIVTGDGMGMALSHGVPLRDMEFVQYHPTGLPGS
NP_418578.1         LVQIRANAVVMATGGAGRVYRYNTNGGIVTGDGMGMALSHGVPLRDMEFVQYHPTGLPGS
                    ************************************************************

WP_001192973.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
NP_418578.1         GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
                    ************************************************************

WP_001192973.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
NP_418578.1         KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
                    ************************************************************

WP_001192973.1      GIETDQNCETRIKGLFAVGECSSVGLHGANRLGSNSLAELVVFGRLAGEQATERAATAGN
NP_418578.1         GIETDQNCETRIKGLFAVGECSSVGLHGANRLGSNSLAELVVFGRLAGEQATERAATAGN
                    ************************************************************

WP_001192973.1      GNEAAIEAQAAGVEQRLKDLVNQDGGENWAKIRDEMGLAMEEGCGIYRTPELMQKTIDKL
NP_418578.1         GNEAAIEAQAAGVEQRLKDLVNQDGGENWAKIRDEMGLAMEEGCGIYRTPELMQKTIDKL
                    ************************************************************

WP_001192973.1      AELQERFKRVRITDTSSVFNTDLLYTIELGHGLNVAECMAHSAMARKESRGAHQRLDEGC
NP_418578.1         AELQERFKRVRITDTSSVFNTDLLYTIELGHGLNVAECMAHSAMARKESRGAHQRLDEGC
                    ************************************************************

WP_001192973.1      TERDDVNFLKHTLAFRDADGTTRLEYSDVKITTLPPAKRVYGGEADAADKAEAANKKEKA
NP_418578.1         TERDDVNFLKHTLAFRDADGTTRLEYSDVKITTLPPAKRVYGGEADAADKAEAANKKEKA
                    ************************************************************

WP_001192973.1      NG
NP_418578.1         NG
                    **

So, in summary, RefSeq will combine multiple identical sequences into a single WP_* multi-species accession. You should therefore expect to find one 100% identical WP_* sequence for your query and multiple, almost identical WP_* entries. And that's precisely what you see here.

Dear terdon, thank you very much for your nice answer. It seems that it's  better for me to do BLAST for this multi-species WP_* accession instead of the specific NP_* accession, isn't it?

@Sara if this is still about finding homologs, then just include whatever you find. The more (relevant) sequences you include, the better your results will be. Now, if the NP and the WP sequences are 100% identical, it makes no difference at all which one of the two you use as a query since they are the same sequence. By the way, feel free to ping me (@terdon)  in the Bioinformatics chat room if you need more clarifications.


text:
What is the distribution of indel sizes in a healthy human genome? of insertion:deletion ratios?
My understanding is that indels are from 1bp to 10Kb, and a healthy genome has ~400K-500K Indels. Surely most of these are small.
What is the distribution of insertion sizes in a healthy human genome? What is the distribution of deletion sizes? What is an average ratio of insertions:deletions?
I have not found a definitive reference.

One of the 2015 papers from the 1000 genomes project has a nice figure (figure 1) showing the size distribution of medium to large sized insertions and deletions:

From another 2015 1000 genomes paper, one can see that the absolute number of smaller indels is much larger, though an exact size range isn't given (as far as I saw). If you really want to know that, just download the most recent 1000 genomes VCF file(s) and compute the exact median size and/or distribution that you want.

Genome-In-A-Bottle (GIAB; version 3.3.2) contains 3.21M SNPs on auto+X chromosomes and 0.51M INDELs in 2.58Gb confident regions. The ins:del ratio is 0.92. On the CHM1-CHM13 pacbio assembly (European ancestry), there are 3.57M auto+X SNPs and 0.58M INDELs in 2.71Gb confident regions [reference] with an ins:del ratio 0.99. These give you an idea of relative short INDELs.
For indels longer than 50bp, I would highly recommend to read this paper by Chaisson et al. This call set is constructed from the whole-genome assembly of the CHM1 genome. It is far more comprehensive and probably more accurate than all the other call sets. The following is Table 1 from the paper:

I have to say that the high ins-to-del ratio is surprising to me, but this is what the data tells us. I have done a similar analysis and reached a similar ratio.
PS: These numbers are essentially from one haplotype. They will be higher for a diploid sample.
EDIT: some numbers in the original post are wrong.

"GIAB ignores hard regions that tend to harbor more indels" -- that's interesting. Do you have a source for that or could you explain more?

@Randoms Sorry. Some of my numbers are wrong. Please see the update. As to your question, GIAB excludes long STRs, while half of indels come from STRs. You can see that "(0.58-0.51)/(2.71-2.58)=0.54" indels per kb outside the GIAB confident region is much higher than "0.51/2.58=0.20" inside its confident region.

Ok thanks. Where are those numbers in that coming from?

@Randoms These VCFs are public. You can count by yourself.

Could you explain where the numbers for the CHM1-CHM13 dataset are coming from? The Syndip paper you linked refers to "0.38 million 2–50bp INDELs". The rep2.37.broad.hc.raw.vcf.gz VCF, by my count, has 936587 indels—493900 are within [2, 50]. The full.37m.vcf.gz VCF has 1145454 indels—556256 in [2, 50].

Intersecting the VCFs with the full.37m.bed.gz, which I guess I should have been doing, yields lower estimates—though still not what quite you're quoting, or what's in the paper strangely.

@Randoms There is a lot to counting: the version of data and manuscript, and the way you take intersection and count. Comments are not the best place for discussions.

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:
How can I find mutations associated with disease in human histone residues?
I would like to look if there are mutations in residues of human histones associated with any disease. For instance, if a mutation in residue K6 (lysine 6) of histone H2A1A is associated with any human disease, (with PUBMED evidence, preferably), but in a systematic way (for all types of histones and all residues).
The purpose of this is to get some candidate residues to study. I have checked the OMIM database but I didn't find anything relevant. It seems that there is no systematic way to connect specific residues to diseases and search for them.
Is there any web database or API available for this purpose?

Hi plat. What kind of database are you looking? One with residues in proteins and diseases? Have you looked at GWAS to find nucleotid sequence and then aa residue? BTW the question in the title is different from the one in the body. Which one want you answered?

I am looking for association of a protein (specifically histones) residue to a disease. So a database of residues in proteins and diseases that I can look at in a systematic way will be fine. However, what you comment about GWAS and check the extrapolated aa is a good idea also. Any recommendation?

Sorry, but I find your question too vague: What have you tried to do? What failed? What other databases has you searched? Why each one failed? What do you lack from those databases?. Please narrow down the question (by [editing]) or ask your advisor how you can approach this.

Would HistoneHits be of use? https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2666297/

The problem with HistoneHits is that it is for yeast not for humans, but thank you anyway!

I have tried OMIM and I didn't find any way of relate histone residues mutations (or nucleotides and its correspondent residues) with human diseases. I don't know if this is because there is no such association or because this is not the correct database to search. Because of that I am asking if there is any database to look at this problem (i.e associate mutations in human histone residues to diseases in human). I think the question is more or less specific enough.

Apologies, I can't read it would seem!

Ensembl contains this information: When you go to the “phenotype” menu item of a given gene, you will see a list of variants (potentially after clicking on “ALL associated variants”) with their associated genomic position.
You can subtract the transcript’s start position from that position to find out the residue (of the translated sequence, and be careful to subtract the lengths of all preceding introns!). Unfortunately this seems to require some manual work.
Alternatively you can query the individual SNPs, which will tell you their position in all transcripts (in this case: Leu 298/68/244):

I found DisGeNET useful for my purpose, a database that associates genes with diseases and, if known, gene variants with diseases. It integrates several sources of information for computing a score of association of gene/variation with the disease ( Piñero et al., 2016, Piñero et al., 2015):

Gene Disease Associations from UniProt, PsyGeNET, ClinVar, Orphanet, the GWAS Catalog, CTD (human data), and Human Phenotype Ontology, RGD, MGD, and CTD (mouse and rat data) GAD, LHGDN and BeFree
Variants Disease Associations from ClinVar, the GWAS Catalog, Uniprot, GAD, and from BEFREE

If available, it provides the Pubmed ID of articles supporting the proposed association and the rs id from the dbSNP for variants associated to disease.
In addition it provides an experimental R package, disgenet2r, that allows to query the database in a programmatic manner and that provides some plot functions to generate representations of the associations founded.
Using this approach, I have been able to found several human histone-disease associations (not my original purpose), together with some concrete histone residues that are mutated in specific diseases.

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


text:
After artificially creating events in a FASTA file, how do I keep track of the old coordinates?
I'm beginning with the reference genome in FASTA format, hg19. I am reading the sequence into a Python dictionary with BioPython:
genome_dictionary = {}
for seq_record.id in SeqIO.parse(input_fasta_file, "fasta"):
    genome_dictionary[seq_record.id] = seq_record

This creates a Python dictionary with keys as chromosome names and values as strings from the FASTA.
My goal is to take the FASTA, and simulate artificial SNPs and InDels along the FASTA to benchmark various algorithms. If I was only creating SNPs in the FASTA, this would be a trivial problem. The reference FASTA has the original coordinates, and the SNPs would also be in a coordinate system which is the same.
However, when I randomly insert an insertion (e.g. an insertion of 12 bp), I move all events into a new coordinate system such that
new coordinate system = old coordinate system + 15

This becomes very difficult to keep track off throughout the entire genome after several interations. I therefore cannot keep track of the new and old coordinate system without running an algorithm and creating a VCF, which may have errors. This defeats the purpose of the simulation.
The output would be either a BED file or VCF which has keep track of the changes.
What method/data structure is used to solve this problem? Would LiftOver work? Maybe CrossMap?
http://genome.sph.umich.edu/wiki/LiftOver
http://crossmap.sourceforge.net/

What are you using to create the changes? Why not just keep track of their location with it?

@DevonRyan These are just string manipulations, e.g. delete 10 characters is a 10 bp deletion. Keeping track isn't trivial----a delete of 10 bp means all bases must be moved over 10 from the original location.

Very informally, I'd just suggest that many bioinformatics file formats just don't really handle changed coordinate systems very well. If you can somehow do the benchmarking without having to manage changed coordinates, you might be better off.

"If you can somehow do the benchmarking without having to manage changed coordinates, you might be better off." This doesn't seem possible

@ShanZhengYang: Err, that is trivial to track, you only care where the changes are in the original file, not where the positions are in the result.

@DevonRyan Perhaps you could explain? I'm randomly inserting SNPs, insertions and deletions in the FASTA. After +3 million events, I need to know how the new coordinates relate to the old coordinates, e.g. a VCF. This sounds quite algorithmically difficult to track to me. How would you do this?

@ShanZhengYang: As you write the the output, track (1) where you are and (2) whether you insert/delete/mutate the base. As you make a modification, write it to a VCF or other file. You don't need to do anything other than that.

@DevonRyan This still doesn't make sense. So, randomly draw all of the positions for the indels/SNPs. Then sort them. Then implement them one by one, keeping track of the changes made?

@DevonRyan "As you make a modification, write it to a VCF or other file." I think these modifications must be in order. Otherwise, it doesn't work. Correct?

@ShanZhengYang: You can define at each base what, if any, kind of variant to add there. Yes, you MUST do this in order. Whether you predefine where the variants will be (modify that from the back then) or do it as you write out the sequence (modify if from the front then) is entirely up to you and makes no real difference.

"with keys as chromosome names and values as strings": I think that the code you show, the values are not simple strings, but SeqRecord objects.

@bli You're correct---I should edit.

@DevonRyan This is making sense now. Thank you for the help. A separate question is if I wanted to simulate a small structural variation between chromosomes, like a chromosomal translocation....in this more complex example, it becomes more difficult to keep track of changes.

@ShanZhengYang: Yeah, don't go crazy with those.

@DevonRyan What do you do with overlapping intervals for events? Here, the order will matter. e.g. insertion and SNP vs. SNP and insertion

Any one location will only ever have one (or zero)  modifications.

@DevonRyan My question may have been unclear. You sort the events. When you do this, it's possible several events will occur at the same (or mre likely, overlapping) genomic coordinates. For example, an insertion from [5,15] and a SNP at 10.

That would be more conveniently described as an insertion at base 5 (10 bases) and an SNP at base 10. Don't think in terms of InDel end positions, just start and widths. Obviously a deletion can't overlap a SNP, so keep track of that.

Note: I see in the documentation that there is a SeqIO.to_dict function that seems to do your data reading in one step

I've written a handful of programs from scratch to simulate mutations and variations in real or simulated sequences.
The trick has always been to sort the variants by genomic coordinate, apply the variant with the largest coordinate first, then apply the variant with the second largest coordinate, all the way down to the variant with the smallest coordinate.
Indels and other variants that affect sequence structure only affect subsequent coordinates. So going in reverse order ensures variants at the beginning of the sequence don't throw off variants at the end of the sequence.

"going in reverse order ensures variants at the beginning of the sequence don't throw off variants at the end of the sequence" This is ingenious, thanks! Remaining question: what would be the most efficient way to create the VCF? @DevonRyan above has recommended the approach of writing each variant into a VCF when you manipulate the sequence. Is this your approach, or would you suggest something different?

That's going to give you the quickest results, because there are existing tools available for creating modified sequences from VCF files.

Separate question: Have you ever tried the above with structural variations? Let's say I wanted to simulate a small structural variation between chromosomes, like a chromosomal translocation....in this more complex example, it becomes more difficult to keep track of changes.

@ShanZhengYang Nope, I've only ever done SNVs, indels, and inversions. Nothing more complex.

If you look for a program which would randomly introduce SNPs + short indels and then would save everything into a VCF file, DWGsim or Mason Variator could be a good choice. Then you can create a corresponding Chain file using bcftools consensus -c and transform various formats between these two coordinate systems using CrossMap.

How do these algorithms generate an accurate VCF after randomly introducing SNPs and InDels?

The SNPs / InDels are generated directly by these programs so the programs know the coordinates, etc.

I don't understand how they keep track of these though. If you randomly insert an InDel, it changes the entire reference sequence.

@ShanZhengYang Mason probably generates an initial read in the reference coordinate. If the initial read overlaps with an indel in VCF, it modifies the read sequence to include the indel. This way, the final read uses the reference coordinate and contains the indel. You don't need to think about another coordinate system at all.


text:
Range overlap python error with genomic regions
I have two files
s3.txt :
1   10  20
1   5   20
2   20  30
2   25  30
1   10  50
2   20  60
1   14  17

s4.txt:
1   10  20
2   20  30

I am trying to match col0 of both the files and get rows that fall between range(inclusive of themselves) 10-20 and 20-30 as seen in s4 file.
file s4 has co ordinates which can be used as reference range (chrom start and end) and s3 has list of co ordinates from an experimental condition what I am trying to achieve is to which co-ordinates from my file s3 fall on or between my reference co ordinates in s4.
code so far:
containing_ranges = []
with open('s4.txt', 'r') as f:
    for line in f:
        fields = line.strip().split('\t')
        containing_ranges.append(fields)

tested_ranges = []
with open('s3.txt', 'r') as f:
    for line in f:
        fields = line.strip().split('\t')
        tested_ranges.append(fields)

for c_range in containing_ranges:
    for t_range in tested_ranges:
        tst = int(t_range[1])
        ten = int(t_range[2])
        cst = int(c_range[1])
        cen = int(c_range[2])
        if  c_range[0] == t_range[0]:
            included = cst >= tst and cen <= ten
            if included == True:
               print t_range

Output with missing row(1   14  17) :
['1', '10', '20']
['1', '5', '20']
['1', '10', '50']
['2', '20', '30']
['2', '20', '60']

Desired output:
1   10  20
2   20  30
2   25  30
1   14  17

Not sure if my logic is wrong and why does it miss 14-17 as it falls between 10-20
[EDIT] using pybedtools

>>> print(s4.intersect(s3, wb=True))
1   10  20  1   10  20
1   10  20  1   5   20
1   10  20  1   10  50
1   14  17  1   14  17
2   20  30  2   20  30
2   25  30  2   25  30
2   20  30  2   20  60

>>> print(s4.intersect(s3, wa=True, wb=True, F=1))
1   10  20  1   10  20
1   10  20  1   14  17
2   20  30  2   20  30
2   20  30  2   25  30

using bedops
bin$ less answer.bed
1       5       20
1       10      20
1       10      50
1       14      17
2       20      30
2       20      60
2       25      30

using @bli code(on python2.7)
('1', 10, 20)
('1', 14, 17)
('2', 20, 30)
('2', 25, 30)
 why can I not see the interval 1 5 20

What result do you get with bedops --element-of?

Please bring your logic/coding questions to Stack Overflow. The relationship of your question to the subject of Bioinformatics is merely coincidental.

@RobertC If OP adds a "bed" tag, this question will immediately look like a bioinformatics question. Also, see the answers. OP is much more likely to get such spot-on answers here. This question could be improved for sure, but it is not off-topic.

Just use bedops, as indicated. Using wrappers to command-line tools is rarely a substitute for learning the tools.

Can you please add more story/context around this question? It has the appearance of a pure programming question (presumably why it was marked as off-topic). It would be nice if you could explain what the different numbers mean, and why you want to do this.

You should use more meaningful names for your variables. It would make the code easier to read, for others but also for you.

It seems you want the ranges from s3.txt ("tested ranges") that are included in a range in s4.txt ("containing ranges"). In this case, I think your error is in the start and end coordinates comparisons. Let's note t_start and t_end the coordinates of the tested range and c_start and c_end the coordinates of the containing range. What you want is c_start <= t_start and t_end <= c_end.

I edited your question to use more meaningful variable names, and I also removed 1    5    20 from the expected output: If I understood correctly, this is not what you want because this is not included in any of the ranges defined in s4.txt

I can't post an answer since your question is "on hold", but here is a (hopefully) corrected version of your code, with minor some coding style improvements and using python3: http://paste.ubuntu.com/24915950/
I hope this helps.

ohh I was away for a bit, so many suggestions thank you all. will edit my post once I go through each of your suggestions

@AlexReynolds added answer

@bli thank you for cleaning up code added answer

You say "which co-ordinates from my file s3 fall on or between my reference co ordinates in s4". If I interpret this correctly, this means that you accept partial overlaps too, not only complete inclusions. Then the desired output should be all the ranges in s3, and not the restricted list I mistakenly corrected.

@novicebioinforesearcher It looks like bedops was able to find your missing interval. If you need to deal with strand labels in the sixth column (per BED specification), you can split a BED file by strand via awk '$6=="+"' in.bed > in.forward.bed and awk '$6=="-"' in.bed > in.reverse.bed, and then run set operations on each of the strand-split files. If you need to reconstruct one file at the end, use bedops -u to do a multiset union of all input BED files.

pyranges answer:
# pip install pyranges
# or
# conda install -c bioconda pyranges

import pyranges as pr
s3 = pr.read_bed("s3.bed")
s4 = pr.read_bed("s4.bed")
s3.intersect(s4, how="containment")

Answer with setup (Edit):
import pandas as pd
from io import StringIO
import pyranges as pr

c1 = """1   10  20
1   5   20
2   20  30
2   25  30
1   10  50
2   20  60
1   14  17"""

c2 = """1   10  20
2   20  30"""

columns = "Chromosome Start End".split()

df1 = pd.read_table(StringIO(c1), sep="\s+", header=None, names=columns)
df2 = pd.read_table(StringIO(c2), sep="\s+", header=None, names=columns)

gr1 = pr.PyRanges(df1)
gr2 = pr.PyRanges(df2)

print(gr1.intersect(gr2, how="containment"))

Result:
+--------------+-----------+-----------+
|   Chromosome |     Start |       End |
|   (category) |   (int32) |   (int32) |
|--------------+-----------+-----------|
|            1 |        10 |        20 |
|            1 |        14 |        17 |
|            2 |        20 |        30 |
|            2 |        25 |        30 |
+--------------+-----------+-----------+
Unstranded PyRanges object has 4 rows and 3 columns from 2 chromosomes.
For printing, the PyRanges was sorted on Chromosome.

You're reinventing bedtools intersect (or bedops), for which there's already a convenient python module:
from pybedtools import BedTool

s3 = BedTool('s3.bed')
s4 = BedTool('s4.bed')

print(s4.intersect(s3, wa=True, wb=True, F=1))

The wb=True is equivalent to -wb with bedtools intersect on the command line. Similarly, F=1 is the same as -F 1.

NotImplementedError: "intersectBed" does not appear to be installed or on the path, so this method is disabled.  Please install a more recent version of BEDTools and re-import to use this method. used pip install pybedtools

You should be able to read and understand how to deal with that error message.

Hi i installed bedtools  made is environment variable i can use it in bash from any directory i still get the same error, is there a different way than this? thank you

You can always install it via conda, though it works for me with the outdated version of bedtools that I have installed (2.25.0).

i was able to get it working thanks, but the result seems to be ambiguous have added it to my original post, 5-20 seems to be missing

Whether you want F=1 or f=1 will depend on your exact needs and the order you do things.

Regarding your python code
If you want the experimental ranges that are entirely contained in one of the reference ranges, you need to have the coordinates in the following order:
cst <= tst < ten <= cen

If what you want are the experimental ranges that overlap one of the reference ranges, you need to have either the start or the end of the experimental range fall within the reference range:
(cst <= tst < cen) or (cst < ten <= cen)

Your code is equivalent to neither possibilities:
cst >= tst and cen <= ten

This is equivalent to:
tst <= cst and cen <= ten

(Or tst <= cst < cen <= ten, since we know that cst < cen, by definition of a bed interval).
With this rewriting, we can more easily see that you are actually selecting the experimental ranges that contain a reference range.
Here is some (python3) code that gives you the results in the other two situations:
#!/usr/bin/env python3

ref_intervals = []
with open("s4.txt", "r") as f:
    for line in f:
        (chr, start, end) = line.strip().split("\t")
        ref_intervals.append((chr, int(start), int(end)))

exp_intervals = []
with open("s3.txt", "r") as f:
    for line in f:
        (chr, start, end) = line.strip().split("\t")
        exp_intervals.append((chr, int(start), int(end)))

contained = []
overlapping = []
for (r_chr, r_start, r_end) in ref_intervals:
    for (e_chr, e_start, e_end) in exp_intervals:
        if e_chr == r_chr:
            if r_start <= e_start < r_end or r_start < e_end <= r_end:
                overlapping.append((e_chr, e_start, e_end))
            if r_start <= e_start < e_end <= r_end:
                contained.append((e_chr, e_start, e_end))

print("overlapping")
for (chr, start, end) in overlapping:
    print(chr, start, end, sep="\t")

print("contained")
for (chr, start, end) in contained:
    print(chr, start, end, sep="\t")

If I run it, I obtain the following results:
$ ./overlap.py
overlapping
1   10  20
1   5   20
1   10  50
1   14  17
2   20  30
2   25  30
2   20  60
contained
1   10  20
1   14  17
2   20  30
2   25  30

It is probably an good exercise to program this, but as the other answers point out, there are efficient tools that would be a better solution in a professional setting.

this makes sense, I guess I need to start thinking and talking like a programmer which would help me write the questions like one.  Sure there are other tools but doing it myself helps me to learn and understand whats exactly going behind the tools.

You could use BEDOPS, instead:
$ sort-bed s3.txt > s3.bed
$ sort-bed s4.txt > s4.bed
$ bedops --element-of 1 s3.bed s4.bed > answer.bed

If you need to run it from within Python, you could use subprocess.check_output():
import subprocess
...
try:
    result = subprocess.check_output("bedops --element-of 1 %s %s > %s" % (set_a_fn, set_b_fn, answer_fn), shell=True)
except subprocess.CalledProcessError as err:
    raise SystemExit("Could not run bedops\n")
# do stuff with 'result'

thank you I would certainly try it , but can you please help me understand what is wrong with my logic/code

You're better off not reinventing the wheel, is my advice. Python is not the right tool for set operations.

If you really want to do this, use debug statements (sys.stderr.write(...)) so you can see what's happening at each stage.

If this is a homework assignment and you are required to write your own solution from scratch, please add a "homework" tag to your question to make this explicit. We're happy to help with homework questions, but the responses will likely be different (and hopefully accommodating).

If it's not a homework assignment, Alex is right: better to use well-tested tools for interval arithmetic.

@DanielStandage yes it is homework will make changes to tag, tried doing it it would not allow me to create a tag

@novicebioinforesearcher I've added the homework tag to your post.

@AlexReynolds does the method suggested account for strand as well?

I'm not sure. I don't see anything in your code that separates intervals by strand. You could use awk ($6=="+") etc. to do this.

@AlexReynolds sorry for late reply i was curious you are right my code does not handle strand, bedops is new for me hence wanted to ask you upfront will install and update you with results


text:
What is a good pipeline for using public domain exomes as controls?
I'm currently attempting association analysis with an extremely small set of patient exomes (n=10), with no control or parental exomes available. Downloading the ExAC VCF of variant sites (http://exac.broadinstitute.org/downloads) or the 1000G integrated call sets (http://ftp.1000genomes.ebi.ac.uk/) and combining this with our pooled patient VCFs has not been successful (I suspect the approach of attempting to merge such large VCFs generated from different pipelines is rather naive).
Looking at the primary literature, I have gathered it should be possible to use these resources to help increase statistical power for our analysis. My question is how do I take these large .vcfs with many samples and successfully merge them to our patient .vcfs, such that the combined VCF can be used downstream to run analysis packages? (PODKAT, PLINK, etc.)

I gather that you want to use the background SNP frequency as a prior for input to your SNP calling algorithm?
I'm not sure of a canned algorithm for doing this, but a quick google shows up some promising links:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3593722/
https://software.broadinstitute.org/gatk/documentation/article.php?id=4723
https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-015-0489-0

Thanks for the response - to clarify, I need additional details on how to correctly assign my priors using these public resources to jointly analyse multiple samples (our patient exomes and the population controls). The first link touches on this, but I'm unsure how to actually implement this approach.
The third link would be useful, except our coverage is fine and instead we lack the sample size they have by several orders of magnitude (gulp)!

Yeah, I don't know of any tools that do this using existing data. What you /can/ do is re-analyze your data in the context of the other data in a single multi-sample call. This seems like a hack TBH, but should work.


text:
How can I find the relevant pathway map from gene-gene or protein-protein interaction list?
I have some difficulties to understand/interpret the pathway map and how a gene-gene interaction list or DNA sequencing can map into pathways.
In addition what's the difference between MARK/ERK pathway and the output of KEGG

You don't build pathway maps from bioinformatics data, you build them from wet-lab experiments.

I agree with Devon about patways because pathways indicate there is a shared known function.
Protein-protein interactions/regulations maps are however another thing where bioinformatics can help out. One of the better exampels is www.string-db.org. Other options are by looking for clusters of correlation/anti correlation in large volumes of expression data.

@0x90 Could you please clarify your question? How do you define a pathway  (The link is to the MAPK/ERK "pathway" in wikipedia)? This is linked with the Kristoffer's observation. And which representation do you want from a pathway (If you want any)?  BTW you seem also interested  in comparing representation of pathways if I understood correctly the last sentence, that would be a different question probably

What do you mean by bioinformatics data?

@Llopis I mean how do you find the affected pathways from a gene-gene or protein-protein interaction list?

That's a very different thing to what your question is asking ("pathway map"). Please [edit] your question and clarify what you are trying to do.

@terdon I am not sure, how to properly rewrite the question.

I can't help, I'm afraid since I don't know what you are trying to do :) A good start would be to show us an example of the data you have and an example of the kind of output you would like from it. Your links to wikipedia seem strange. What are we supposed to take from that? Also, gene-gene interactions (gene regulatory networks) are not the same as protein-protein interactions (which are closer to what rection pathways are). So it is really hard to understand what you're asking.

See a project like Plant reactome, where they infer the pathways of all sequenced plant genomes using orthology to well annotated species.
http://plantreactome.gramene.org/

As you alluded to in your question, KEGG does provide a curated set of enzyme and metabolite information. This information can be parsed and used to create a network that you can analyze to look at how gene products could be working together. Additional steps could include mapping transcriptome data to the map.
See this paper as an example of how this can be done.
I hope this helps. :)


text:
what is the difference between local identity and global identity for homolog finding?
We believe that if after running blast, the global identity between a resulting sequence and our query is at least 30%, we can say that two sequences are homologs. what is the difference between local identity and global identity and how can we calculate them?
the following file is the global alignment result of one of my PSI-BLAST hits against query:
Aligned_sequences: 2
# 1: NP_418577.1
# 2: WP_036312822.1
# Matrix: EBLOSUM62
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 264
# Identity:      85/264 (32.2%)
# Similarity:   125/264 (47.3%)
# Gaps:          29/264 (11.0%)
# Score: 376.5
#
#
#=======================================

NP_418577.1        1 -----------MAEMKNLKIE-------VVRYNPEVDTAPHSAFYEVPYD     32
                                ..|..:..|:       :.|::||||..|....|:|...
WP_036312822.      1 MTIVDSGAPADTQEANDSGIQSYLVTFIIRRFDPEVDAEPRWVDYDVEMY     50

NP_418577.1       33 ATTSLLDALGYIKDNLAPDLSYRWSCRMAICGSCGMMVNNVPKLACKTFL     82
                     .|..:||||..||.::...||:|.||...||||..|.:|...:|||||.:
WP_036312822.     51 PTDRVLDALHRIKWDVDGTLSFRRSCAHGICGSDAMRINGRNRLACKTLI    100

NP_418577.1       83 R--DYTDGMKVEALANFPIERDLVVDMTHFIESLEAIKPYIIGNSRTADQ    130
                     :  |.:..:.|||:...|:|:||:|||..|.||...::|::...|.....
WP_036312822.    101 KDLDISKPIYVEAIKGLPLEKDLIVDMDPFFESFRDVQPFLQPKSAPEPG    150

NP_418577.1      131 GTNIQTPAQMAKYHQFSGCINCGLCYAACPQFGLNPEFIGPAAITLAHRY    180
                     ....|:....|.|...:.||.|..|.::||.|..:.::.|||||..|||:
WP_036312822.    151 KERFQSIKDRAVYDDTTKCILCAACTSSCPVFWTDGQYFGPAAIVNAHRF    200

NP_418577.1      181 NEDSRDHGKKERMAQLNSQNGVWSCTFVGYCSEVCPKHVDPAAAIQQGKV    230
                     ..||||.....|:..||.:.|||.|.....|:|.||:.::...||.:.|.
WP_036312822.    201 IFDSRDDAADVRLDILNDKEGVWRCRTTFNCTEACPRGIEITKAIAEVKQ    250

NP_418577.1      231 ESSKDFLIATLKPR    244
                     ...:.
WP_036312822.    251 AVLRG---------    255

#---------------------------------------
#---------------------------------------
the following file is the local alignment of that hit against query:
Aligned_sequences: 2
# 1: NP_418577.1
# 2: WP_036312822.1
# Matrix: EBLOSUM62
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 219
# Identity:      83/219 (37.9%)
# Similarity:   119/219 (54.3%)
# Gaps:           2/219 ( 0.9%)
# Score: 390.5
#
#
#=======================================

NP_418577.1       13 RYNPEVDTAPHSAFYEVPYDATTSLLDALGYIKDNLAPDLSYRWSCRMAI     62
                     |::||||..|....|:|....|..:||||..||.::...||:|.||...|
WP_036312822.     31 RFDPEVDAEPRWVDYDVEMYPTDRVLDALHRIKWDVDGTLSFRRSCAHGI     80

NP_418577.1       63 CGSCGMMVNNVPKLACKTFLR--DYTDGMKVEALANFPIERDLVVDMTHF    110
                     |||..|.:|...:|||||.::  |.:..:.|||:...|:|:||:|||..|
WP_036312822.     81 CGSDAMRINGRNRLACKTLIKDLDISKPIYVEAIKGLPLEKDLIVDMDPF    130

NP_418577.1      111 IESLEAIKPYIIGNSRTADQGTNIQTPAQMAKYHQFSGCINCGLCYAACP    160
                     .||...::|::...|.........|:....|.|...:.||.|..|.::||
WP_036312822.    131 FESFRDVQPFLQPKSAPEPGKERFQSIKDRAVYDDTTKCILCAACTSSCP    180

NP_418577.1      161 QFGLNPEFIGPAAITLAHRYNEDSRDHGKKERMAQLNSQNGVWSCTFVGY    210
                     .|..:.::.|||||..|||:..||||.....|:..||.:.|||.|.....
WP_036312822.    181 VFWTDGQYFGPAAIVNAHRFIFDSRDDAADVRLDILNDKEGVWRCRTTFN    230

NP_418577.1      211 CSEVCPKHVDPAAAIQQGK    229
                     |:|.||:.::...||.:.|
WP_036312822.    231 CTEACPRGIEITKAIAEVK    249

#---------------------------------------
#---------------------------------------

BLAST (Basic Local Alignment Search Tool) doesn't give you global alignments, so you may need to reconsider what you're doing/asking.
Furthermore, "homolog" has a very specific definition. Homologs are genes that are both descended from the same ancestral sequence. It doesn't, as many often mistake, mean 'similarity'. 30% sequence identity, does not imply homology by itself. You need more information.

Take a look here (http://www.ebi.ac.uk/Tools/psa/) if you want to find tools for local and global alignment.

no it's not my problem. Maybe I couldn't describe it well. I have run PSI-BLAST for my sequence and I 've got some sequences that I believe they are homologs of my sequence, but I'm not sure if they are really homologs or not and which one of these sequences are real homologs. For making sure about this homology, I must check local and global identity and I had done it. Maybe a part of a sequence with mine, has a good identity, but I can't say these two sequences are homologs. They must have globally good identity

Perhaps I'm wrong but "They must have globally good identity" still isn't really enough to say that these sequences are homologs. You need ancestry information for the organisms you're interested in.
If you really want to know orthology/homology, I'd suggest using a dedicated orthology finding program like OrthAgogue, OrthoMCL, or better yet ROARY if your sequences are bacterial.

No they don't. It all depends on how closely related the species are. Please [edit] your question and add more information. There is no magix sequence simlarity threshold that defines homology.

how much the species should be close?

There is no magic formula. For closely related species I would expect a sequence identity of far more than 30% for a true homolog but for distant species, I could also expect far lower. As I told you in my answer to your previous question, it is far more important to look at the conservation of specific residues than to look at the overall sequence similarity value.

Please don't post images of text. They just make the question slower to load, they cannot be edited and cannot be searched. Instead, just copy/paste the text directly into your question and use the formatting tools to make it appear like code. Also, the alignment is useless if we don't i) know the species being aligned and ii) know the e-value of the alignment.

If you want a purely alignment based approach to begin to infer homology (note I say infer, as we stated in the comments, this isn't sufficient to say 2 sequences are homologs), I would use a Multiple Sequence Alignment (MSA) instead of a pairwise like you've done so far, and include as many similar proteins as you can from BLAST.
The choice of algorithm for the MSA probably won't matter too much, as you will be looking for conserved regions anyway.
Once you have your MSA, you can easily spot conserved domains. If anything is known about this protein structurally, it's quite trivial to work out whether a conserved region is a functional domain (active site, that kind of thing).
With this information, you would be somewhat justified in calling these proteins structural homologs at least if they share domains, but, to reiterate, without doing phylogenetic analyses, you cannot call them true homologs.

Defining homology isn't as simple as an identity cutoff, however, you're right that homology seems likely given a certain identity. This is where e-value comes into play. See:
https://en.wikipedia.org/wiki/Expected_value
To answer the question, global-global requires that the alignment account for every base in the query and the target sequence. local-global requires a the full length of the query on a part of the target. local-local allows only the highest scoring part of the alignment between the query and target to be reported. i.e. one functional domain against another. these hits are often called HSPs.


text:
How to represent trans-spliced genes in GTF?
For example, see this gene (nad1) in ENA:
http://www.ebi.ac.uk/ena/data/view/ABI60879
If you look at the XML for that gene you see the following:
join(
             DQ984518.1: 324706 .. 325091 ,
  complement(DQ984518.1:  24417 ..  24498),
  complement(DQ984518.1:  22828 ..  23019),
             DQ984518.1:   3484 ..   3542 ,
  complement(DQ984518.1: 153702 .. 153960)
)

Which shows 5 exons joined out of phase and out of order. Is there a valid GTF representation of this?
How to dump a 'non-canonically spliced' gene into GTF? i.e. what's the recommendation?
Cross-posted on biostars

Maybe you could explain better what is the particularity of splicing for this gene. It is not easy to see in the link you give.

Yes, that's much clearer. I don't have an answer, though

To my knowledge there's no defined way to deal with that in GTF. GFF3 handles trans-splicing (you'll have to scroll down to "trans-spliced transcript") by giving an individual transcript multiple parents (e.g., ID=some_transspliced_gene;Parent=gene1,gene2). You could use the same methodology with GTF files, but just note that it'll break most downstream programs.

To be fair, annotating this in spec-compliant GFF3 would break the vast majority of GFF3 parsers, too. :-)


text:
How to dump genes from GenBank in GFF3 format?
This question has also been asked on BioStars
If I look at this record in GenBank I see about 6k genes:
https://www.ncbi.nlm.nih.gov/nuccore/CM000760?report=gbwithparts
I'd really like to be able to dump those genes in GFF3 format, but I'm guessing I can't do that in a single step, so...

How to dump all genes from a GenBank record in 'whatever' format?
How to convert 'whatever' format to GFF3?

Turns out you can just grab the GFF3 from the NCBI's FTP site!
ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/003/195/GCA_000003195.3_Sorghum_bicolor_NCBIv3/
See:

https://www.ncbi.nlm.nih.gov/genome/doc/ftpfaq/#files

Or "Access the data" on the right here:

https://www.ncbi.nlm.nih.gov/assembly/GCF_000003195.3

That's useful, but it isn't actually answering the question you asked. Could you please either change the question so it is asking where you can get the data (instead of how to convert GenBank/anything to GFF) or delete either your answer or the question? I know this is kind of weird, but these are strictly Q&A sites, so answers should be answering the question asked.

When accessing all of the annotated genes for a reference genome, downloading a GFF3 file directly from the Genbank or RefSeq FTP sites is definitely the way to go.
But when submitting a search query into NCBI Entrez, there is no easy way to download the results in GFF3 format. Probably the best route is to download the results in ASN.1 format, and then use annotwriter from the NCBI C++ toolkit to convert to GFF3. But getting the latter to compile and run properly isn't trivial in my experience.


text:
Find gene at position from gff or gbk file
I have a VCF file with SNPs from a bacterial genome and want to find if the SNPs are located inside genes, is there some CLI-tool where you can pass a VCF file and a gff or gbk file and it returns the name of the genes?

Please show us an example of your inputs. What format are your SNPs in? Are they in a proper VCF file? Are they in some other format one per line? Many per line? Do you want a command line solution? an online graphical tool? Please [edit] your question and clarify.

Via BEDOPS:
$ gff2bed < annotations.gff > annotations.bed
$ vcf2bed < snps.vcf > snps.bed
$ bedmap --echo --echo-map-id-uniq snps.bed annotations.bed > answer.bed

This can be reduced to a one-liner if you're using bash:
$ bedmap --echo --echo-map-id-uniq <(vcf2bed < snps.vcf) <(gff2bed < annotations.gff) > answer.bed


text:
Correlating gene expression with qualitative variables
I have a gene expression dataset that I want to investigate. Particularly, I would like to understand whether there is any correlation between each gene's expression and some quantitative or qualtitative data (say, correlation between gene 'XPTO' , body mass index, and race).
One possible way to test this would be through logistic regression, but is this a good approach or are there caveats that I should know about using such a statistic?
My question is the following: which methods would you advise to measure such correlations, and why?
(this post was crossposted on Biostars)

Would a linear regression not make more sense as expression is not a binary outcome?

Logistic regression would generally be a bad choice for non-binary outcomes. In such cases, linear regression (or a GLM more generally) still works fine. You can already do that in the standard R packages for RNAseq (DESeq2, edgeR, and limma), where the fold-change is in whatever units you're measuring your quantitative trait in.


text:
Best distance parameter for estimating physical interaction between residues in a PDB file
We can calculate the distance between residues in a PDB file regarding different parameters like closest atoms, alpha carbon, beta carbon, centroid and etc. Which one of these parameters are better to show physical interaction between residues in a PDB file?

Which kind of physical interaction are you looking for: an hydrogen bound, a covalent bond, a ionic bond...?

There is no different for me which kind of physical interaction. The least distance and the farthest no different just physical

The CASP experiment uses an 8 Angstrom distance between C-beta atoms (C-alpha for glycine) to determine residue-residue contacts.

Use a distance cutoff of 12Å between Cᵦ atoms.

If this is part of your previous answer, consider editing it, otherwise it would be nice to provide bit more explanation.

Please could you clarify what types of interactions can happen at such distances, as I am currently unaware of interaction being possible at such a long distance?

A popular (and I would say, respected and trusted) website is PDBsum https://www.ebi.ac.uk/pdbsum (which also has a Wikipedia article about it: https://en.wikipedia.org/wiki/PDBsum)
They measure protein-protein (prot-prot) contacts as any Nitrogen, Carbon, or Oxygen element atoms (N, C, O) within 4.0 Angstroms (Å) when measuring the 3d euclidean distance.
If your atomic contacts are not close enough, the amino acids may not actually interact at all, even though they theoretically could.  Since in PPIs you generally want to use reliable ground-truth information, it is better to use cautious parameters.  Often, in PPIs, the answer to "does it interact?" is not "yes" or "no" but rather "yes" or "don't know / unknown".

You should read my paper, where I did a comprehensive parameter sweep across different thresholds and definitions and benchmarked the resulting empirical potentials on a decoy ranking challenge:
The relevant phrase from the abstract is "compared using 90 different definitions of residue-residue contact". However, it boils down (as always) to the question, "what do you want to do".
Cite me if you use my paper for your research.

Most papers give 5 or 6 angstrom distance for physical interactions. I don't know this distance is between which atoms exactly

People use what other people use, which is why I did a 'parameter sweep'. I've seen other work using a strict 'information maximization' criteria, which was nice.


text:
qPCR: Why is fold change and standard deviation calculated after transformation?
I am analyzing data from a quantitative polymerase chain reaction (qPCR) using R. After cleaning the raw data, it looks something like this:
> dput(x)
structure(list(Reporter = c("FAM", "FAM", "FAM", "FAM", "FAM",
"FAM", "FAM", "FAM", "FAM", "FAM", "FAM", "FAM", "VIC", "VIC",
"VIC", "VIC", "VIC", "VIC", "VIC", "VIC", "VIC", "VIC", "VIC",
"VIC"), Number = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L,
12L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L), A = c(22.19,
22.24, 22.5, 22.54, 22.6, 22.59, 23.07, 23.46, 22.43, 22.74,
24.09, 23.91, 24.52, 25.03, 25.25, 25.82, 25.13, 24.71, 25.34,
25.85, 25.25, 26.15, 25.81, 25.29), B = c(21.72, 21.78, 22.86,
22.73, 19.88, 20.07, 21.06, 21.06, 20.96, 21.11, 19.46, 19.43,
24.75, 24.69, 25.64, 25.19, 23.76, 23.69, 24.35, 25.05, 24.1,
23.81, 22.81, 23.13), C = c(21.37, 21.56, 20.07, 20.01, 21.17,
21.08, 20.54, 20.36, 33, NA, NA, NA, 23.91, 24.31, 23.61, 23.88,
24.33, 24.31, 23.37, 23.53, 33, NA, NA, NA), E = c(26.26, 27.33,
25.93, 26.56, 25.76, 23.03, 24.72, 25.27, 24.43, 24.31, 26.98,
23.33, 24.04, 25.02, 25.1, 25.1, 24.68, 25.48, 25.87, 26.22,
25.35, 25.36, 25.11, 25.98), F = c(25.81, 26.9, 25.58, 26.61,
25.06, 21.85, 23.59, 24.04, 23.19, 23.19, 25.17, 20.8, 24.12,
24.26, 25.32, 25.25, 24, 23.78, 24.7, 24.48, 23.52, 23.87, 23.05,
23.05), G = c(26.12, 27.02, 24.08, 25.15, 25.99, 23.18, 24.2,
24.05, 33, NA, NA, NA, 23.47, 23.45, 23.7, 23.74, 24.46, 24.19,
23.56, 23.53, 33, NA, NA, NA)), .Names = c("Reporter", "Number",
"A", "B", "C", "E", "F", "G"), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, 24L))

I think that each well measures two genes: target which is 12 different genes (color FAM), and the reference or housekeeping gene, GAPDH (color VIC). Also, control is triplicate A, B, C (3 x 12 wells), and treatment is E, F, G (3 x 12 wells).
I have created a function ddCt_ for analyzing the data with the ddCt algorithm (Livak & Schmittgen, 2001). It takes one argument x which is a data.frame of the form exemplified above.
ddCt_ <- function(x) {
  # Subset x by control/treatment and target/reference
  # Then, calculate Ct averages for each triplicate
  TC <- x %>%
    select(Reporter, Number, A, B, C) %>%
    filter(Reporter == "FAM") %>%
    rowwise() %>%
    mutate(Ct_Avg = mean(c(A, B, C), na.rm = TRUE))

  RC <- x %>%
    select(Reporter, Number, A, B, C) %>%
    filter(Reporter == "VIC") %>%
    rowwise() %>%
    mutate(Ct_Avg = mean(c(A, B, C), na.rm = TRUE))

  TT <- x %>%
    select(Reporter, Number, E, F, G) %>%
    filter(Reporter == "FAM") %>%
    rowwise() %>%
    mutate(Ct_Avg = mean(c(E, F, G), na.rm = TRUE))

  RT <- x %>%
    select(Reporter, Number, E, F, G) %>%
    filter(Reporter == "VIC") %>%
    rowwise() %>%
    mutate(Ct_Avg = mean(c(E, F, G), na.rm = TRUE))

  # Normalize Ct of the target gene to the Ct of the reference gene
  dCt_control   <- TC$Ct_Avg - RC$Ct_Avg
  dCt_treatment <- TT$Ct_Avg - RT$Ct_Avg

  # Normalize dCt of the treatment group to the dCt of the control group
  ddCt <- dCt_treatment - dCt_control

  # Calculate fold change
  fc <- 2^(-ddCt)

  # Calculate avg and sd
  dCt_control_avg   <- mean(dCt_control)
  dCt_control_sd    <- sd(dCt_control)
  dCt_treatment_avg <- mean(dCt_treatment)
  dCt_treatment_sd  <- sd(dCt_treatment)

  # Create output
  df <- data_frame(
    Sample = 1:12, "dCt Control" = dCt_control, "dCt Treatment" = dCt_treatment,
    "ddCt" = ddCt, "Fold Change" = fc, "dCt Control Avg" = dCt_control_avg,
    "dCt Control SD" = dCt_control_sd, "dCt Treatment Avg" = dCt_treatment_avg,
    "dCt Treatment SD" = dCt_treatment_sd
    ) %>% round(2)
  write.csv(x = df, file = "result.csv")
  saveRDS(object = df, file = "result.rds")
  df
}

However, I am not a biochemist or molecular biologist, so I am unsure of the basic concepts involved and the overall approach. So, my questions are:

Is the function correct?
Why is it important to calculate fold change and standard deviation after normalization?

May I mention ReadqPCR and NormqPCR Bioconductor packages designed exactly for this:
https://www.bioconductor.org/packages/release/bioc/html/ReadqPCR.html
https://www.bioconductor.org/packages/release/bioc/html/NormqPCR.html

The functions look correct, but calculate a few by hand and ensure they match. One thing I should note is that the subtraction of the Ct values usually happens before an average is made, since generally both are in the same well. In other word, subtract the reference Ct from the gene of interest Ct from the same well and then average the dCts. This order is important since one usually does a multiplexed qPCR.
You have to normalize to a reference gene to control for how much cDNA was used, since that will alter the Ct values. If you calculated the fold-changes without normalization then they could be purely due to using more/less cDNA in the reaction (i.e., the output would be meaningless).

Thanks for your reply. I don't think calculation by hand will help since it is not the calculations that I am unsure of, but the overall algorithm: the order of calculations and the type of arrangements and so on. And, could you perhaps elaborate on your reply to my second question.

@Samuel: I've expanded my reply and noted that one normally averaged the dCts, rather than averaging the Cts.

I think that each well measures two genes: target which is 12 different genes (color FAM), and the reference or housekeeping gene, GAPDH (color VIC). Also, control is triplicate A, B, C (3 x 12 wells), and treatment is E, F, G (3 x 12 wells). Does this information give you any further clues? I added this information in the description as well.

@Samuel: Yes, that's quite helpful. It's indeed a standard multiplexed prep. Do make sure to subtract the Cts from the same wells then.


text:
Intersection of two genomic ranges to keep metadata
I am trying to find intersection of two genomic ranges (gr1 and gr2) and keep metadata from one of them
gr1
chrI [1, 100] * | 0.1
chrI [101, 200] * | 0.2

gr2
chrI [50, 150] + |

intersect(gr1, gr2) will remove metadata, where
subsetByOverlaps(gr1, gr2) return:
chrI [1, 100] * | 0.1
chrI [101, 200] * | 0.2

where metadata is saved, but intersection is bigger than real (the same as in initial gr1 start/end)
How can I get an intersection range and keep metadata?
to get the range:
chrI [50, 100] * | 0.1
chrI [101, 150] * | 0.2

Presuming you don't have cases of multiple ovelaps
> library(GenomicRanges)
> gr1 = GRanges(c("chrI", "chrI"), IRanges(start=c(1, 101), width=c(100, 100)), mcols=data.frame(foo=c(0.1, 0.2)))
> gr2 = GRanges(c("chrI"), IRanges(start=c(50), width=c(101)))
> o = findOverlaps(gr1, gr2)
> grl1 = split(gr1[queryHits(o)], 1:length(o)) # You can't mendoapply on a GRanges object
> grl2 = split(gr2[subjectHits(o)], 1:length(o))
> foo = function(x, y) {
+     rv = x
+     start(rv) = max(start(x), start(y))
+     end(rv) = min(end(x), end(y))
+     return(rv)
+ }
> unlist(mendoapply(foo, grl1, y=grl2))

The output is:
GRanges object with 2 ranges and 1 metadata column:
    seqnames     ranges strand | mcols.foo
       <Rle>  <IRanges>  <Rle> | <numeric>
  1     chrI [ 50, 100]      * |       0.1
  2     chrI [101, 150]      * |       0.2
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

It'd be more memory efficient to iterate over o, but this gives you the idea.
In cases of multiple overlaps, one could use intersect(), then findOverlaps() with that and then continue on with the procedure demonstrated above (or iterate over the overlaps).


text:
How are PDO and PDX used in computational and predicative models for tumour biology?
There are wet methods: Patient-derived models: Patient Derived Xenograft (PDX) and Patient Derived Organoids (PDO) to reflect tumor biology.
Are there any databases or computational tools that use the outcomes from PDO/PDX experiments to create a predictive computational model for cancer or other diseases?

What do you mean by outcomes of PDO/PDX experiments? RNA-seq of the tumor? Clinical variables along time for PDX experiments? Interactions between cell lines in the PDO experiments? What do you want to predict: mortality, growth, some clinical variable value, metastasis...? AFAIK they are used to learn about the tumors, not to predict anything.

@Llopis I think that all what you have just listed could use to build a predictive model.

But what do you want to predict ? And more important, would that be a good model for what?

@Llopis by prediction I mean to build a general model to predict if a patient has cancer and which medicine would be the most suitable for him.

You can't predict if a patient has cancer from a cancer model! If you want to predict which medicine for a given patient's tumor would be the most suitable for him you need to test (in silico or not) all the medicines you think they have a change. (Also you would have already removed the tumor or part of it for both methods, so presumably the patient would no longer have it given the current techniques). In my opinion, this question is too broad, if you have a concrete example or a paper doing something similar it could be answerable but as it is, is too broad and lacks many details.

@Llopis saying too broad is a nice way of saying not feasible?

It is feasible, but this question is not answerable.

@Llopis https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4323263/

Let us continue this discussion in chat.

@Llopis I am experiencing some difficulties to login. Please accept my apologies. I am looking for something similar to pan-cancer, other cancer types and Alzheimer's. moreover as you mentioned that's only PDX.

It is Ok. Do your on due diligence and search in the literature for it, when (if) you find yourself stuck ask for the concrete problem you are in. But avoid asking a broad general question.

@Llopis thanks, I still don't think it's too broad, but appreciate your feedback and will try to update it.

this question asks 2 quite unrelated questions, it would be best to separate these into 2 questions.

@nsheff do you mean the first and the secondly? They are closely related.

@0x90 There's no good mechanism to split a question other than to edit it and post a new one. Having said that, your first question is appropriate for the biology SE and only the second one for here.

Potential pitfall!
I'm not sure about predictive models, but you need to be aware of a potential pitfall in blindly aligning PDX or PDO based sequencing data without first removing contaminating host organism reads, as otherwise these will lead to a lot of false positive variants caused by miss alignment.  In my experience even a small mount of host material can lead to a ten fold increase in called variants due to miss aligned host reads looking like true variants.  I'd recommend using Xenome, source.
Note Xenome was designed for DNA sequencing it's not splice-aware, so for RNA-Seq there is no optimal solution other than removing reads which align to the host organisms genome.  Although the issue here is that in conjunction with a splice-aware read aligner synteny between the host and grafted genomes might create some interesting problems.  However for RNA-Seq provided FACS or similar confirms low-levels of host contamination I expect levels of expression will not be badly affected.  Although this really needs investigating.
Finally and rather annoyingly Xenome produces none-standard FASTQ so you'll need to fix it's output with:
awk '{if (NR % 4 == 1) print \"@\"$0; else if (NR % 4 == 3) print \"+\"$0; else print $0 }' as reported here on seqanswers.

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


text:
Publicly available, free, complete database for antibiotics names and classes?
This is a tough one I think: is there a publicly available, up-to-date, free, complete database for antibiotics names and classes?
I am specifically looking for information like, e.g., cefpirome (is a) cephalosporin.
I've looked in a couple of different places, but nothing really comes close. E.g.:

CARD comes closest, but is not free (companies need to pay a license) and not complete/up-to-date (e.g., cefpirome is missing)
@Karel mentioned KEGG in his answer, but this also needs a license (companies for sure, academics if they download data via FTP).
Novel Antiobiotics Database has no class ... and stops in 2003
The Antimicrobial Index is web only, and terms of use apparently prohibit any kind of web-scraping.
A suggestion by @Pierre on Twitter was to try Antimicrobial stubs on Wikipedia. Good idea, but the class is not stored in the ChemBox and must be free-text parsed. Or writing a script which evaluates the Wikipedia category tags.
I drew blank at different government agencies (EMA, FDA, EFSA, NIH, etc.). The closest I found was from the WHO a list of critically important antibiotics. But not complete, and in a PDF.

Any pointer appreciated.

The FDA's NDC may be a good resource. It is freely available and was updated: 6/22/2017. You'll have to do some grepping to get the info you need but I think it is mostly all there. Cefpirome is not in the db but many cephalasporins are. There is also an unfinished drug list which I did not check.

Good one, didn't know it. But unfortunately error prone. Some entries contain multiple ingredients, e.g., "Erythromycin and Benzoyl Peroxide", which fall into different categories and cannot be re-assigned unambiguously with the table data.

Can you define what classes you want? For example, are efflux pumps a class? Classification schemes are not universally; I suggested CARD on Twitter because it's a full ontology rather than just 2 columns

Well, e.g., I'd like to have the association that cefmenoxime is a cephalosporin (bonus points for 3rd, 4th or 5th generation, but not a must). Twitter didn't show me a notification of your response until I specifically looked for it, sorry. But as I wrote in examples: CARD is not free, companies need a license.

ChEBI
It should be possible to get this information from the ChEBI database, see the exported tables. You could download the ontologies (in OWL / OBO), parse them using some ad-hoc parser or using a dedicated library (e.g., Pronto), build a directed acyclic graph based on the is_a relationship, and extract the antibiotic subtree.
KEGG
You could also use KEGG: Kyoto Encyclopedia of Genes and Genomes, but a subscription is probably required. The exported drug table should be easy to parse (e.g., in Python).
For a quick verification of the completeness, you can download the table, filter out everything except antibacterial agents and check the number of records:
curl "ftp://ftp.genome.jp/pub/kegg/medicus/drug/drug" > keg_all.txt

cat keg_all.txt \
   | perl -pe 's/\n/BREAK/g' \
   | perl -pe 's@///BREAK@///\n@g' \
   | grep 'CLASS       Antibacterial' \
   | perl -pe 's/BREAK/\n/g' \
   > keg_antibacterial.txt

cat keg_antibacterial.txt | grep "///" | wc -l

There are 646 records. However, some of them are almost identical (e.g., only 1 atom difference).
If you add | grep -i "antibiotic" \ to the pipeline, then you get only 445 records.

Excellent resource, but does not pass the "freely available" filter: KEGG requires a yearly license. I will update my question to include that.

Would ChEBI satisfy your definition of "freely available"?

ChEBI would, indeed! And it contains approximately what I need. OBO or OWL are not necessarily formats from which it is easy to extract the wanted information as one needs to build and traverse graphs to piece together required data, but it is doable. Can you please amend your answer to include ChEBI before KEGG? I'll accept the answer then.


text:
How to maximize fastq parsing with FastqGeneralIterator (Bio.SeqIO.QualityIO)
I'm contributing to a python-based project that uses Biopython to analyze fastq files. It currently uses SeqIO.parse, which populates various structures with all of the fastq information (including converting quality scores). There is apparently a faster (lighter-weight) parser called FastqGeneralIterator that doesn't populate all of these items.
I'd like to use FastqGeneralIterator, but be able to perform some of the functions SeqIO.parse offers on the fly (not for every sequence). For example, I'd like to convert base quality to PHRED scores for specific sequences, but I don't see that function available.
Is there an easy way to take a tuple output by FastqGeneralIterator and convert it into a proper Bio.SeqIO SeqRecord?

Not sure where is your problem exactly. Just in case it might be of interest to you, you can find some examples of python code parsing fastq (including FastqGeneralIterator, and faster alternatives) in the last part of this answer: https://bioinformatics.stackexchange.com/a/380/292

You may also be interested in this: http://biopython.org/DIST/docs/api/Bio.SeqIO.QualityIO.FastqPhredWriter-class.html (There's a FastqPhredWriter in Biopython)

I had not come across FastqGeneralIterator before, but I will start using it for some of my work!
One answer is to replicate the code in the FastqPhredIterator, which does exactly what you are looking to do.
However, My personal preference would be to use the biopython infrastructure as much as possible, especially when it comes to the tricky aspects of record parsing.  If the records you want extra info on are uncommon, then using FastqPhredIterator for phred quality extraction might work for you:
import io
from Bio.SeqIO.QualityIO import FastqGeneralIterator, FastqPhredIterator

with open("input.fastq", "rU") as handle:
    for title, sequence, quality in FastqGeneralIterator(handle):
        # do something with title, sequence, and quality
        # Set special_case to True when you want more info on the record
        special_case = False
        if special_case:
            record_stream = io.StringIO("\n".join([title, sequence, "+", quality]))
            record = next(FastqPhredIterator(record_stream))
            # do something special with record

What we are doing here is combining the record strings into one and dumping them into a StringIO object, which implements the file handle duck-type. We then hand the in-memory handle to FastqPhredIterator to parse the FASTQ record again. Make sure to set special_case to True only when you want extra info on a record, as the special_case branch makes at least two extra copies of the FASTQ record.

Very useful, thank you! I was hoping there were dedicated methods for these tasks.


text:
PDB format: remark number for free text
I would like to add a text to PDB files that I'm processing with my tool, rna-pdb-tools. Someone points that the way I'm using it right now it's not correct (https://github.com/mmagnus/rna-pdb-tools/issues/48).
I use HEADER right now which is not correct.
HEADER Generated with rna-pdb-tools
HEADER ver 37c5b4e-dirty
HEADER https://github.com/mmagnus/rna-pdb-tools
HEADER Mon Oct 10 22:18:19 2016
ATOM      1  P     C A   1     -19.687  -3.296  65.469  1.00  0.00           P

Do you happen to know which remark number to use? (http://www.wwpdb.org/documentation/file-format-content/format33/remarks1.html)?

From this document, it looks like REMARK 250 is the way to go.

Great! I decided to use your answer https://github.com/mmagnus/rna-pdb-tools/commit/b11d33eba8e1311e89d98ba70e2eb087db878982

Why? This remark is reserved specifically for details of the experimental method (e.g.  powder diffraction) that the model is based on.

EDIT:
Actually, after having a more thorough look into the documentation; REMARK 3 might be the more relevant to what you need. Other tools are reported in that section.

REMARK 3 presents information on refinement program(s) used and related statistics. For non-diffraction studies, REMARK 3 is used to describe any refinement done, but its format is mostly free text.

More on Remark 3

From the same source as you reference yourself v 3.3, I would argue that
REMARK 0 would be more accurate than REMARK 250.

REMARK 0 (updated), Re-refinement notice
  REMARK 0 identifies entries in which a re-refinement has been performed using the data from an existing entry. This remark also describes the PDB code and the journal records for the original data set.

While when you consider the definition of REMARK 250
-- which follows directly:

REMARK 205 (specific to Fiber diffraction experiment),
REMARK 201(and REMARK 215/217 all specific to NMR experiment),
REMARK 230 (neutron diffraction study),
REMARK 240(electron crystallography study) and
REMARK 245(and REMARK 247: specific to EM study)

REMARK 250, Other Type of Experiment Details
REMARKs specific to other kinds of studies, not listed above.
  REMARK 250 is mandatory if other than X-ray, NMR, neutron, or electron study. The format of the date in this remark is DD-MMM-YY. DD is the day of the month (a number 01 through 31), MMM is the English 3-letter abbreviation for the month, and YY is the year.

and then REMARK 265 is also about the crystallography experiment.
Hence, I know that personally, I usually assume that whatever is annotated in REMARK 250is relevant to the first acquisition experiment to me.
Another reason I am more inclined to think REMARK 0is more accurate, it is due that the definition of REMARK 0 was updated while REMARK 6-99 are are no longer for use of free text annotation.
BTW, I am not sure how Pdb_extract would do with the remarks when converting the PDB file to CIF files. It would probably worth perusing  PDB to PDBx/mmCIF Data Item Correspondences page.

Use any REMARK number not reserved for something else.
The REMARK field was initially for general remarks. In the format description from 1992 only the numbers 1 to 3 were reserved for specific things, with the note:
Other general commentary is given in higher numbered REMARKs.
Then, in the next years, other remark numbers got reserved (and sometimes retired) for various specific annotations. The content of remarks was getting more and more formalized, to make it easier to parse by computers. The specifications from 1996 and 1998 recommended:

Non-standard remark annotations, or those with no clearly-defined
  topic or assigned remark number, appear with remark number 6 or
  greater, but less than remark number 100.

I'd still follow this advice.
It is generally understood that when parsing PDB files one may have custom REMARK records.
I've seen also programs that output remarks without numbers - and this also doesn't cause problems with any reasonable PDB reader.
On the other hand, using remarks with numbers reserved for very specific things
(such as remark 250 in the most upvoted answer) can cause confusion.
The current spec (from 2012) is a precise description of files served by the wwPDB. It doesn't attempt to define a format for other uses, such as passing data between programs. In the section about REMARKs 6-99 it has the last bullet point saying: "REMARKs 6-99 are no longer for use of free text annotation", because these remarks were removed from the old PDB entries. But it doesn't prevent from using it outside of the wwPDB.


text:
Which quality score encoding does PacBio use?
Do you know which quality score encoding PacBio uses now? I know some of their file formats have changed in the past year or two, but I haven't found much on their quality score encoding.
The most recent answer I found is from 2012, where one user states it's "phred style", but doesn't state whether it's 33 or 64: DeNovo assembly using pacBio data

PacBio does use PHRED 33, but it turns out the question may be irrelevant for the newer PacBio Sequel Sequencer, because it reports all base qualities as PHRED 0 (ASCII !). The RS-II reports PHRED 33 quality scores.
The Sequel provides data in the form of subreads, which are the circular consensus sequences (CCS) from a single zero-mode waveguide (ZMW). Essentially, a single circularized DNA molecule enters a ZMW well, and is sequenced as many times as possible. (Sometimes multiple molecules get into a single well, but that's not the normal case.) The number of times this single molecule can be sequenced depends on several variables, including molecule length. The consensus sequence for this single molecule is then reported out as a subread.
As I started playing with the quality scores in my data, I realized the quality score for every base is !, which equates to a PHRED quality of 0. I contacted someone at PacBio who informed me that they learned the base quality values were not any more useful than the bases themselves for generating a high-quality circular consensus sequence. Thus, rather than spend the computational effort to calculate the scores, they skipped it, and reported ! for every base. They also report an overall read quality of 0.8, but this is simply a placeholder. They are working to make this information more readily available.


text:
NaN values after ComBat analysis on TCGA COAD RNAseq
I have FPKM-UQ data from COAD-TCGA.
I generated an expression set of this data using:
> edata = log2(data + 1)
> edata[1:2,1:2]
      X01240896.3f3f.4bf9.9799.55c87bfacf36
1                                  8.540967
10                                13.528968
100                               16.296422
1000                              13.658199
10000                             15.143788
      X01ad5016.f691.4bca.82a0.910429d8d25b
1                                  9.886537
10                                16.719682
100                               17.212312
1000                              10.317842
10000                             13.166767

Row name is entrezid, and column names are case_ids.
The batch data provided by TCGA which indicates a total batch num = 42
> omics[1:5,]
                                      snf    boo order batch
X01240896.3f3f.4bf9.9799.55c87bfacf36   3 192807     1    22
X01ad5016.f691.4bca.82a0.910429d8d25b   2  82213     2     2
X01f493d4.229d.47a6.baa8.32a342c65d01   1 A08907     3    29
X022f39e9.57ee.4b2b.8b3a.8929e3d69a37   2 A16S13     4    32
X02f9668c.71e6.485f.88b1.b37dc8bdd2ab   3 102113     5     7
> batch = omics$batch
> batch
  [1] 22  2 29 32  7  1  8  8  7 41 26 14 26 26  1  5 22  5 15 18 28 12 12  8 19
 [26] 13 38  5 27  7 41  1 22 26 28 22 33 14 13  8 10 13 11 17  5 38  5  7  7 31
 [51] 22 42 11  8  7 17 15 17 12 17  2 19  2 12 18 14 22 11 21 15 17 17 33 15 19
 [76] 15  6  1 10 11  2 28 36 17  9 15 12 17 15 10  6  8 13 25 11 15  1 13 24  8
[101] 13 15 17  8  2 22 17 11 13 37 19 38 13 19  5 16  5  5 22 35 35 19 13 17 22
[126]  7  8  8 41 39 14  6 19 18 38 22 14 30 19 24 11 14 13 16  7 13  8 22 12  8
[151] 22 17  2 40 13 35 15  5 24 19 41 22 22 19 17  8 13 33 10 17 29 12  9 11 24
[176] 27  5 14  7 25 13 35  9 18 16 13 31 25 18 10 35 18 14 33  5 19 17  8 11 11
[201] 42  7 10  8 39 19  8  7 11 19 31 35 28 36 13 37 10 38  3 18 13 34 39 13  3
[226] 40 31 17 19 26  7 19 12 22 19 17 12 17 13 12 28 33 37 17 35 18 40 13 20 35
[251] 17  6 15 14 12 27  5 41 31 41 18 21 19 17 40  8 41  2  5 15  5 15 15  5 39
[276]  1 22 33 13 17  7 11 15  3  7 38 15  7 14 14 22 31 14 18  3 19  7  9  6  5
[301] 22  8 11  1 13 10 19  8 14 15  6 32  1  2 19  4  5 13 18 42 11 17  8 36 19
[326] 15 22 28 15  8 11  8 15 38 19 39  7 11 42 23  5  1 22 17 35 13 40 25  7 41
[351] 38  5  8 19 31  6 38  8 36  4 15 15  7 12 22 38 29 35 15  2 10  2 19 15 25
[376]  5 39 30 19 18 15 10 10 33  1 19 32 19 19 15 20 15  7 22 35 22 30 14 12 16
[401] 21 13 15 18  7  5  6 31 12  7 15 32 33 18 14 15 15 26 31  1 14 19 17 18  1
[426]  5 10 32

Unfortunately, when I run the ComBat analysis it generates all NANs, and no values:
> modcombat = model.matrix(~1, data=omics)
> combat_edata = ComBat(dat=edata, batch=batch, mod=modcombat, par.prior=TRUE, prior.plots=FALSE)
> combat_edata[1:5,1:2]
      X01240896.3f3f.4bf9.9799.55c87bfacf36
1                                       NaN
10                                      NaN
100                                     NaN
1000                                    NaN
10000                                   NaN
      X01ad5016.f691.4bca.82a0.910429d8d25b
1                                       NaN
10                                      NaN
100                                     NaN
1000                                    NaN
10000                                   NaN

What happens if you make batch a factor?

So when I normalize the data first by mean, then perform the log2 conversion, the batch correction works.

Stipulating batch as factor generates the same NAN result. I'm thinking that maybe it wasn't close enough to a normal distribution without the mean normalization, so it wasn't able to calculate the correction?

Hi Martin, thanks for your question and welcome to Bioinformatics Stack Exchange. I have changed your title to be more specific and match the question that you have asked. Having specific questions makes it more likely that people will be able to answer your problem (or be interested in answering it), and some people won't click through if they see a question that has a vague title.

It's just interesting because the data is supposedly already normalized between samples using FPKM-UQ. So I'm not clear why the additional mean normalization is necessary. In particular, I need to do something like: data <- standardNormalization(data) to scale the mean at 0 w/ stdev @ 1, then log2(data + 1).

I think I might know what the problem is.
I plotted the log2 corrected values, and the distribution may not be normal because of the 0 values present in RNASeq data. If I understand correctly, this data is actually missing, whereas with my current analysis it is viewed by R as an actual value of 0. So either the missing data needs to be imputed, or the data needs to be scaled and normalized to mean of 0, then subsequently log2 transformed.
> data <- standardNormalization(data)
> data <- log2(data + 1)
> combat_data = ComBat(dat=data, batch=batch, mod=modcombat, par.prior=TRUE, prior.plots=FALSE)
> combat_data[1:5,1:2]
          X01240896.3f3f.4bf9.9799.55c87bfacf36
100507661                           -0.13857428
57103                               -0.06700506
22838                                0.16295096
55567                               -0.13579991
6147                                 1.64990772
          X01ad5016.f691.4bca.82a0.910429d8d25b
100507661                           -0.13789534
57103                               -0.06210232
22838                                0.01546224
55567                               -0.13443289
6147                                 1.48473072


text:
how to set database other than nr for remote blast+ search
I am attempting to run a BLAST search remotely using BLAST+.
I can get search to work correctly at the command line with the following commands:
blastp -query proteins.fasta -remote -db nr -out proteins_nr.txt -outfmt 6 -evalue 1e-30

However, I would like to us the remote database titled "Microbial proteins from nr" which is what would be used for the microbial blast search accessing BLAST from the NCBI site.
I'm struggling to find the correct code to access this database instead of just "nr".
Is there a list with the codes of different databases?
Thanks in advance

I believe you're looking for env_nr? It's listed as such, under Metagenomic proteins in the blastp webpage. It appears that the word within brackets should be supplied alongside the -db parameter. A quick test with a dummy amino acid fasta file does turn up a result to a valid NCBI protein accession.

Edit: I followed up on this with a little more digging around, and it appears the database you're looking for is called Microbial_proteins. I ended up quote-searching Microbial proteins from nr on Google and exactly 2 hits turn up. The first hit, I think, is what the OP posted a screenshot of in a comment to this answer. The second is a list of databases along with keywords from a Spain-based bioinformatics firm: http://data.biobam.com/ncbi_blast_dbs_protein.pdf. I retried the dummy example with the Microbial_proteins database, and it seems to work. Again, I am not sure why I can't easily find this on an NCBI resource. Also, I am curious as to how different this database is from env_nr.

@flatley176 thank you for response. That answers what I can specify with -db. Unfortunately I'm not looking to blast against env_nr. I would like to blast against the website database using blast+ and specify this database:

is there a way to specify this for a -remote blast using blast+
thank you.

I do not know how Microbial proteins from nr differs from the env_nr database. Page 3 of the BLAST manual (ftp://ftp.ncbi.nlm.nih.gov/pub/factsheets/HowTo_BLASTGuide.pdf) lists the available databases - the only reference to the word "microbial" is to the 16s rRNA. That said, I dug around a little more, and it appears there IS a solution - the database is called Microbial_proteins. Will update my answer.


text:
How do you write a .gz fastq file with Biopython?
How do you write a .gz (or .bgz) fastq file using Biopython?
I'd rather avoid a separate system call.
The typical way to write an ASCII .fastq is done as follows:
for record in SeqIO.parse(fasta, "fasta"):
    SeqIO.write(record, fastq, "fastq")

The record is a SeqRecord object, fastq is the file handle, and "fastq" is the requested file format. The file format may be fastq, fasta, etc., but I do not see an option for .gz.
Here is the SeqIO API.

I'm not sure I'm doing it the best way, but here is an example where I read a compressed gzip fastq file and write the records in block gzip fastq:
from Bio import SeqIO, bgzf
# Used to convert the fastq stream into a file handle
from io import StringIO
from gzip import open as gzopen

records = SeqIO.parse(
    # There is actually simpler (thanks @peterjc)
    # StringIO(gzopen("random_10.fastq.gz").read().decode("utf-8")),
    gzopen("random_10.fastq.gz", "rt"),
    format="fastq")

with bgzf.BgzfWriter("test.fastq.bgz", "wb") as outgz:
    SeqIO.write(sequences=records, handle=outgz, format="fastq")

You shouldn't need the StringIO bit for parsing, just make sure you open the gzip handle in text mode under Python3:
records = SeqIO.parse(gzopen("random_10.fastq.gz", "rt"), "fastq")

@peterjc I didn't know about the "text mode", thanks. That explains the error message I got:ValueError: Is this handle in binary mode not text mode?

Here is how I was able to export to a standard gzipped FASTQ file using Biopython. Basically, instead of using SeqIO.write, I directly called the .format method of the SeqRecord object. The example code below imports a gzipped FASTQ file, removes reads that do not contain a G in positions 7, 8, and 9, and writes the results to a gzipped FASTQ file.
from Bio import SeqIO
import gzip

path_in = "path/to/in.fastq.gz"
path_out = "path/to/out.fastq.gz"
handle_in = gzip.open(path_in, "rt")
handle_out = gzip.open(path_out, "wt")

fq = SeqIO.parse(handle_in, "fastq")
for read in fq:
    # Only export reads that have a G in positions 7, 8,
    # and 9
    if read.seq[6:9] == "GGG":
        handle_out.write(read.format("fastq"))

handle_in.close()
handle_out.close()

Using the code example from the original question that converts a FASTA file to a FASTQ file, it would look something like the following:
fastq = gzip.open("path/to/out.fastq.gz", "wt")
for record in SeqIO.parse(fasta, "fasta"):
    fastq.write(record.format("fastq"))
fastq.close()

Caveat: This only works for Python 3. Python 3 differs from Python 2 in how it imports data from binary files (source). In fact, ideally you should be using Python 3.4 or greater since Biopython has compatibility issues with Python 3.3 (source).

As far as I know you can not write compressed fastq on the fly. You will have to write the entire fastq and then compress it.

Hi Aaron Berlin, thanks for your answer and welcome to Bioinformatics stack exchange. Your answer is a bit short; do you have any source you can provide about not being able to write compressed fastq with Biopython, or an alternative python hook that could be used to work around the inability for it to support writing out a compressed file?


text:
How to obtain .bed file with coordinates of all genes
I want to get a .bed file with the genes' names and canonical coordinates, also I would like to have coordinates of exons, too. I can get the list from UCSC, however, if I choose UCSC Genes - knownCanonical, I can not extract coordinates of exons. If I use other options - I am getting coordinates of as many transcriptional isoforms as were detected while I need only one canonical form.
How can I get such BED file?

So if we have a case with two exons that respecitively have position 1-10 and 5-15 - what coordinates do you want in that case?

Via Gencode and BEDOPS convert2bed:
$ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "gene"' - \
    | convert2bed -i gff - \
    > genes.bed

You can modify the awk statement to get exons, by replacing gene with exon.
To get HGNC symbol names in the ID field, you can add the --attribute-key="gene_name" option to v2.4.40 or later of convert2bed. This slight modification extracts the gene_name attribute from the annotation record and puts it in the fourth (ID) column:
$ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "gene"' - \
    | convert2bed -i gff --attribute-key="gene_name" - \
    > genes.bed

BEDOPS: https://github.com/bedops/bedops
This is based off an answer I wrote on Biostars, which includes a Perl script for generating a BED file of introns from gene and exon annotations: https://www.biostars.org/p/124515/#124522

No need to specify the input “file” - for gunzip.

For those who want to know the last release number (now it should be 27) you may check here: http://www.gencodegenes.org/releases/current.html . So the first URL should be ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_27/gencode.v27.annotation.gff3.gz

The latest link should be ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz, as Sanger no longer participates in Gencode.

On could also load the gff3 file as a dataframe, and save the columns used in the bed file as a tab separated file without headers and index.

You can indeed get this information from the UCSC Table Browser. Select knownGene as your primary table, make a filter, add knownCanonical to as a linked table to filter on, then in the free-form query section add "1" without the quotes. Then click submit and select Bed output, where you can choose "exons plus" as an output option. This will lead to output like the following:
chr1    17368   17436   uc031tla.1_exon_0_0_chr1_17369_r    0   -
chr1    29553   30039   uc057aty.1_exon_0_0_chr1_29554_f    0   +
chr1    30563   30667   uc057aty.1_exon_1_0_chr1_30564_f    0   +
chr1    30975   31097   uc057aty.1_exon_2_0_chr1_30976_f    0   +
chr1    30365   30503   uc031tlb.1_exon_0_0_chr1_30366_f    0   +
chr1    34553   35174   uc001aak.4_exon_0_0_chr1_34554_r    0   -
chr1    35276   35481   uc001aak.4_exon_1_0_chr1_35277_r    0   -
chr1    35720   36081   uc001aak.4_exon_2_0_chr1_35721_r    0   -

For more step-by-step information, please check the answer from the mailing list archives here (in particular the Exon Method 1 section):
https://groups.google.com/a/soe.ucsc.edu/d/msg/genome/BJ-6DlaZNCY/grgGIpuJAwAJ
Also note that UCSC provides several mailing lists for support, if you have futher questions please send them there:

General questions: genome@soe.ucsc.edu
Questions involving private
data: genome-www@soe.ucsc.edu
Questions involving mirror sites:
genome-mirror@ose.ucsc.edu

You can download a list of transcript annotations as a flat file from UCSC:
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/refGene.txt.gz
It is not exactly a BED file, but it does contain information about the known transcripts for this assembly (hg19, in this case):
585 NR_046018   chr1    +   11873   14409   14409   14409   3   11873,12612,13220,  12227,12721,14409,  0   DDX11L1 unk unk -1,-1,-1,
585 NR_024540   chr1    -   14361   29370   29370   29370   11  14361,14969,15795,16606,16857,17232,17605,17914,18267,24737,29320,  14829,15038,15947,16765,17055,17368,17742,18061,18366,24891,29370,  0   WASH7P  unk unk -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,

The gene name can be found in column 13. The exon start and end positions are in columns 10 and 11, respectively. More information can be found in the other columns:

bin
name
chrom
strand
txStart
txEnd
cdsStart
cdsEnd
exonCount
exonStarts
exonEnds
score
name2
cdsStartStat
cdsEndStat
exonFrames

(The column information was obtained from the sql file here)