text:
Download program for ΔG for TM insertion
I am trying to calculate the ΔG of TM insertion for many protein sequences.
Has anyone come across the source code or a programmatic version for the von Heijne ΔG for TM insertion? The web server is at http://dgpred.cbr.su.se/index.php?p=instructions
I can't find a link in the papers or on the site.

Ask the authors otherwise you'll require a "digger", (web spider). I'll share my "digger" code for another delta G website but I have not used this specific site. I think it is common for web servers not to distribute the "backend" for delta G calculations. If you get the standalone please post back. I'm about to start a very similar project (not on TM, but stability of a single antigen) and coding to explore the delta G parameter space is of a lot of interest

@MichaelG. Contacting the authors is the next step. I just want to check that I'm not missing their repo. A digger/scraper isn't ideal, unfortunately. The tool I'm building needs to be a bit of futureproof.

The application is stored in github (ElofssonLab/dgpred) with a link to the website.
To find it I searched for: dgpred download github


text:
Filter VCF for "CAF=[0.9,0.1]"
Related to this other question
but in this case the CAF is defined as "An ordered, comma delimited list of allele frequencies...".
I tried using vcffilterjdk and the command suggested by Pierre
java -jar dist/vcffilterjdk.jar -e 'variant.variant.getAttributeAsStringList("CAF","1.0").stream().skip(1L).map(A->Double.parseDouble(A)).anyMatch(A->A<0.01);' input.vcf

but predictably returned an error
java.lang.NumberFormatException: For input string: "0.1]"

I don't know javascript so how should I change the command to make it work?
EDIT:
HEADER
##INFO=<ID=CAF,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on 1000Genomes, starting with the reference allele followed by alternate alleles as ordered in the ALT column. Where a 1000Genomes alternate allele is not in the dbSNPs alternate allele set, the allele is added to the ALT column.  The minor allele is the second largest value in the list, and was previuosly reported in VCF as the GMAF.  This is the GMAF reported on the RefSNP and EntrezSNP pages and VariationReporter">

LINE
chr1    51479   rs116400033     T       A       .       PASS    ASP;CAF=[0.8926,0.1074];COMMON=1;G5;G5A;KGPROD;KGPhase1;OTHERKG;RS=116400033;RSPOS=51479;SAO=0;SSR=0;VC=SNV;VLD;VP=0x050000000005170016000100;WGT=1;dbSNPBuildID=132

Can you show us the line from your VCF file that causes the error? And also the header line where CAF is defined?

Printing a bracketed list of numerical values is not typical for VCF—I'd guess it's not strictly allowed. A list should be printed only with the commas to separate values. Perhaps the solution is as simple as removing the square brackets.
cat input.vcf | perl -ne 's/CAF=\[([^\]]+)\]/CAF=$1/; print' > fixed.vcf

Replace cat with gunzip -c if your VCF file is gzip-compressed.
The Perl regular expression is

looking for the string CAF=[...]
grabbing the characters in ...
replacing the entire match with CAF=...

Here's a detailed breakdown of the Perl regular expression.

s/.../.../; — substitution syntax; replace the contents of the left block with the contents of the right block
[^\]]+ — one or more non-] characters; this is the trickiest bit; see this tutorial for more info on the [] syntax for character classes and the special ^ symbol for character class negation
(...) — the parentheses capture the matched pattern into the $1 variable
CAF=$1 — replaces everything matched in the first block with CAF= followed by the captured characters

That's probably a better (and easier) solution than changing the javascript code. Thanks for the suggestion.

source: Stack Exchange

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

sort: votes


text:
Restriction enzymes with multiple recognition sequences
In a bioinformatics programming course I am currently teaching (using Python, and approached from a mostly CS rather than a biology background) I assigned a project which involves finding sites where restriction enzymes bind to a DNA sequence. The project is heavily influenced by Chapter 9 of the now classic text "Beginning Perl for Bioinformatics" by James Tisdall. I am using the REBASE database, in particular REBASE Format #9. In this format, each restriction enzyme has a recognition sequence which looks like e.g. C^YCGRG where ^ gives information about the exact cut-point and upper-case letters such as Y and R other than ACGT are ambiguity codes. The point of the project is to turn these into regular expressions where cleavage information is discarded and the ambiguity codes are expanded into character classes, and then use these regular expressions to search for matches.
All fairly straightforward, but one very observant student noticed that some enzymes had multiple sequences, for example M.BssHIII has a total of 5 sequences: ACGCGT,CCGCGG,RGCGCY,RCCGGY,GCGCGC. It is easy to see that these correspond to a total of 11 sequences (after the Rs and Ys are expanded), far fewer than the 144 sequences which would match VSSSSB (which is the least ambiguous single sequence capable of matching the above).
My question is, what is going on here? Am I safe in simply replacing , by | (for "or") in the regular expression, or is there something more subtle going on? Only 8 of the nearly 5000 enzymes have such multiple sequences, and neither the Tisdall book nor even the description of the format (both of which seem to assume single sequences) shed any light on the situation.

It's programming exercise. Just use "|" for regex. It would be also biologically correct.

I would honestly just remove the non-homologous cutting site enzymes from the database. Firstly, these are biology students and require a simple picture to grasp the concept of coding via regex pattern matching. Secondly, in a wet-lab nobody uses these enzymes, I was not even aware the "8" existed until your post. In a teaching lab the best motto is "keep it simple".

I was tempted to remove these 8 from the database, and might end up doing so, but I was genuinely impressed with the student who noticed it and told her that I would look into it. Most of the students are actually CS students taking it as an elective, but I have a couple of biology majors or students that have some major-minor combo between CS and biology or chemistry.

I guess the cat is out of the bag. The difference is this is teaching rather than research. In my experience I would create a subproject to avoid distraction, and offer to all students is interested. Thus if you introduce the | and use it to filter the "anomalous 8 restriction enzymes" into a dictionary, panda, or np, to be analysed later if there is interest. Its interesting a core-bioinformatics text has a "bug".

At first I thought that these enzymes postdated the Tisdall book (which was published in 2001) but a check of the references for these shows that some of them were studied as early as the 1980s. I'll have to check out the book again, it is possible that there was something I missed. It is also possible that in an earlier edition of REBASE these were either missing (since when Tisdall wrote there was only around 1000 enzymes included) or had ? for their recognition sequence. I'll see if it is possible to download earlier versions of it. It could be interesting to see how it changed over time.

Keep in mind 1980s the way cleavage sites were determined different than today. Currently, it's easier to find degeneracy in enzyme cleavage site by NGS. You could also email Dr. Rich Roberts (Nobel Laureate)  and Dana Mecelis. I am sure they will get back to you.

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


text:
Hits in Primer-BLAST not found with programmatic blastn query
Summary:
I'm writing an automated tool for designing Sanger sequencing primers (other tools don't do exactly what I want) that will be incorporated in a production diagnostic laboratory for on-demand NGS variant verification. I'm stuck on a problem wherein off-target BLAST hits found by NCBI's Primer-BLAST aren't showing up in my programmatic BLAST query of the same genome assembly.
Background:
After each iteration of design with primer3, I'm using blastn -task blastn to search a local hg38 (patch 12) database for off-target hits. I have logic that does round-robin pairwise comparisons (i.e. compares every hit to every other hit) to determine if each hit pair could result in an amplicon (on same chromosome, opposite-stranded, perfect match at 3'-ends, plus-strand hit is upstream of minus-strand hit, are separated by <= 1000 bp). If my logic finds a viable off-target amplicon, it triggers a total re-design with different design parameters. If no off-targets are found, those primers are finalized, so it's important that I find any and all off-targets.
To make sure that code for viable off-target amplicons is working correctly (that it's triggering re-design or reporting a design failure when appropriate), I'm using NCBI's browser-based Primer-BLAST to verify the absence of viable off-targets, which also uses genome build hg38 patch 12.
The problem:
Primer-BLAST is finding hits that I can't find with my code. Even when I look at every hit returned from the blastn query, I can't find hits that Primer-BLAST is finding. I would like to know if there is some parameter that I am missing, or parameters I could use in addition to those I'm already using, to produce the hits I am getting from Primer-BLAST.
I believe my blast database an my query fasta file are built correctly, because I am getting logical output from the query below.
Programmatic BLAST Query Details:
Query file (query.fasta):
>for_primer
GGATCTGCCTCCTTGCC
>rev_primer
CCCAACAGGGCCCTC

BLASTn Query Statement:
blastn -db blastn -db /path/to/hg38.fasta -evalue 100 -num_threads 8 -outfmt 6 std gaps nident -dust no -word_size 7 -query /path/to/query.fasta

The BLAST database used hg38.fa.gz (which I believe to be patch 12) from here. That was unzipped and processed with makeblastdb.
Primer-BLAST Results:
The Primer-BLAST results are here. To re-create, use the template DNA below, search database = "Genomes for selected organisms...", organism = Homo sapiens, max target size = 1000. The template DNA is:
 GAGGAGTCAGGCCTGATCTGAAGGCAGTGGGGTACCCTCTAGGGTGTCAGCAGGGCTGTGATGTCACCAGTCTTGGCTGCAGTGTGGATTAGGTAGGCGGAGGCTGGGAGGCAGTGGGAAGTGGCCGCAGGGACCGTGTGGTCAGTCAGGGTGGTGCCATTTCCTGTGTGCACACTCGGGGGCTCCCGCCACACCAGGGTGCAGACCACCCCGCCGCTCTCCCGCCAGGCAGCCTGCAGGGTGAAGGACTGACATGGCGGGAGGAGTGGCCTTCTCTGTTGGGGCGCCTGTGGCCCCCCAGCTCTTTGATAGGCAGAGCTGCAGGATCTGCCTCCTTGCCCCTGTGAGCCCAGCAGTGATTCTTGTGTCCTCCCCAGATCCATGGCTGCGGTGAGAGGGCTCGCGCCCATGCCCGGCAGAGAATCAGCCGGGAAGGGGTCCTCTATGCCGGCAGCGGGACCAAGAACGGATCCCTGGACCCAGCCAAGAGGGCCCAGCTGCAGAGGAGGCTGGATAAGAAGCTGAGTGAGCTCAGCAACCAGAGGACCAGCCGGAGGAAGGAGAGGGGGACGTGACCGGCCGCATCCTTGCACGCCCCGCGGAGCTCTCTCCATGGTAGCCCAGGGCGCTGGCAGACCATGCTCCGCCTCCACCAGGGCCACAGAGGAGCGGAGGGGCCTATGGGGGAGGAGCGGAGGGCCCTGTTGGGGAAGGTTGGGTTTTTGGACCCCAGGGATAAGCTTTTCCGATGTCACAATGTGGAGGAAAGCACCTGGGGGACAACAGTGCTCGTGCAGGTGGGGCTTGGGAAATGCACGTCCCTTCCCCTTACTCCCCGCCAAAACCCACATCCCAGCCTCTGGATCCTGGGGAAGGTTCCAGTCCCTGGAGAATACCCAGGGCCTCAAACTTGAAGTCACTCCTCCAATGTCTGGGACTTGCCAGCTCAGCCCGTTAGGATGAGGGTGCTGAGAGGAAACAGGAAACAAGACTGCGAATGGCG

The blastn output from the query above does not contain either of these two Primer-BLAST hits for 'rev_primer' that comprise a potential amplicon from chromosome 11:
product length = 110
Features flanking this product:
   12759 bp at 5' side: neurexin-2-beta isoform X3
   735 bp at 3' side: RAS guanyl-releasing protein 2 isoform X5

Reverse primer  1         CCCAACAGGGCCCTC  15
Template        64726566  AGT............  64726552

Reverse primer  1         CCCAACAGGGCCCTC  15
Template        64726457  ....G..........  64726471

Why aren't either of these two hits returned in my blastn query? Thanks in advance for helping me with this.
Update:
I noticed in IGV that off-target sequence above is lower-case in the hg38 reference genome. I've tried passing -soft_masking false to blastn, but got the same result.
I thought maybe the problem was due to my using hg38 from UCSC instead of that from NCBI, so I downloaded hg19 from this page and made the blast database from that file. Same result.
I then unmasked the hg38.fna file by converting all lowercase nt to uppercase, as recommended in the NCBI genome assemblies README.txt, created a BLAST database from that and ran my code again using the unmasked database. Same result, still am not detecting the off-target hit above.

I've found the problem, and it's way less interesting than I initially thought. Essentially, my "eValue" was way too low. Instead of -evalue 100, I'm now using -evalue 5000 (Primer-BLAST actually uses 30,000).
For some primers, -evalue 5000 gives ~20,000 results. It's beyond the scope of this question, but I'm now filtering on the query/subject Tm for each hit, perfect match at the query 3'-end, and number of query/subject mismatches. My pairwise
hit comparisons are now detecting plausible off-target amplicons that Primer-BLAST is detecting.
Essentially, I was relying on NCBI's hit ranking system to determine what hits were relevant, but I needed cast a wider net with blastn, and then implement a rational ranking scheme of my own.


text:
Definition of "seed" in sequence alignment
I would like to know what is meant by "seed" for various sequence aligners. How is it important?

The seed is the subset of a read used in the first step of an alignment. Many aligners work by a seed-and-extend model, wherein they first find all regions matching the "seed" and then extend the alignment around that allowing mistmatches and indels until it either gives up (and therefore uses a different seed) or finds a sufficiently good alignment.

Thanks devon for the explanation. I have a slightly similar question as to why finding a seed for a very short exons is difficult? Considering a 3 nt exon, is it because seed  such small has a probability of mapping every 4^3=64 bases and virtually impossible to look for in the genome? Your views on this

Usually seeds are a few times that length, since you want few enough hits to limit your search space quickly.

Devon's answer gives a good, concise definition. But it's also helpful to consider why seed-and-extend is used and what benefits it provides.
Finding approximate string matches requires operations that are expensive computationally. On the other hand, finding exact string matches can be performed using much less expensive operations—that is, it can be done more quickly with less memory. So quickly finding the location of all seed matches and then applying the approximate matching ONLY AT THOSE LOCATIONS (the seed-and-extend strategy) makes sure that less time is wasted on chasing spurious similarities.
So how long should a seed be? Riffing off your comment, a seed 3bp in length would be entirely ineffective since that same 3bp sequence is going to occur ALL OVER THE GENOME IN RANDOM PLACES: genes, pseudogenes, transposable elements, regulatory sequences, you name it. On the other hand, a seed 51bp in length would also be ineffective since almost every seed would contain at least one sequencing error or variant/polymorphism, thus preventing a perfect match.
In practice, mappers designed for aligning 100-300bp Illumina reads typically use seeds in the high teens or low 20s.


text:
I get the following error when using the GSEA software 'Enrichment Map Visualization' tool using the Cytoscape 3.7.1 application
I am trying to visualize a GSEA analysis I conducted using the Broad Institute's GSEA Enrichment Map Visualization, however I get the following error and am not sure what the issue is: any help would be greatly appreciated.
INFO - successfully connected to cytoscape rest
INFO - successfully connected to cytoscape rest
INFO - Found enrichment map command
INFO - Get URL:http://localhost:1234/v1/commands/enrichmentmap/gseabuild?edbdir=C%3A%5CUsers%5CHatim%5Cgsea_home%5Coutput%5Cmar15%5Cmy_analysis.Gsea.1552686798457%5Cedb&pvalue=0.005&qvalue=0.1&overlap=0.5&similaritymetric=OVERLAP&combinedconstant=0.5&expressionfile=C%3A%5CUsers%5CHatim%5CDesktop%5CJunior%5Csem+2%5CKarras+Lab%5CGSEA%5CKarrasEXP.gct.txt
DEBUG - CookieSpec selected: default
DEBUG - Auth cache not set in the context
DEBUG - Connection request: [route: {}->http://localhost:1234][total kept alive: 0; route allocated: 0 of 2; total allocated: 0 of 20]
DEBUG - Connection leased: [id: 4][route: {}->http://localhost:1234][total kept alive: 0; route allocated: 1 of 2; total allocated: 1 of 20]
DEBUG - Opening connection {}->http://localhost:1234
DEBUG - Connecting to localhost/127.0.0.1:1234
DEBUG - Connection established 127.0.0.1:51968<->127.0.0.1:1234
DEBUG - Executing request GET /v1/commands/enrichmentmap/gseabuild?edbdir=C%3A%5CUsers%5CHatim%5Cgsea_home%5Coutput%5Cmar15%5Cmy_analysis.Gsea.1552686798457%5Cedb&pvalue=0.005&qvalue=0.1&overlap=0.5&similaritymetric=OVERLAP&combinedconstant=0.5&expressionfile=C%3A%5CUsers%5CHatim%5CDesktop%5CJunior%5Csem+2%5CKarras+Lab%5CGSEA%5CKarrasEXP.gct.txt HTTP/1.1
DEBUG - Target auth state: UNCHALLENGED
DEBUG - Proxy auth state: UNCHALLENGED
DEBUG - http-outgoing-4 >> GET /v1/commands/enrichmentmap/gseabuild?edbdir=C%3A%5CUsers%5CHatim%5Cgsea_home%5Coutput%5Cmar15%5Cmy_analysis.Gsea.1552686798457%5Cedb&pvalue=0.005&qvalue=0.1&overlap=0.5&similaritymetric=OVERLAP&combinedconstant=0.5&expressionfile=C%3A%5CUsers%5CHatim%5CDesktop%5CJunior%5Csem+2%5CKarras+Lab%5CGSEA%5CKarrasEXP.gct.txt HTTP/1.1
DEBUG - http-outgoing-4 >> Host: localhost:1234
DEBUG - http-outgoing-4 >> Connection: Keep-Alive
DEBUG - http-outgoing-4 >> User-Agent: Apache-HttpClient/4.5.2 (Java/1.8.0_201)
DEBUG - http-outgoing-4 >> Accept-Encoding: gzip,deflate
DEBUG - http-outgoing-4 >> "GET /v1/commands/enrichmentmap/gseabuild?edbdir=C%3A%5CUsers%5CHatim%5Cgsea_home%5Coutput%5Cmar15%5Cmy_analysis.Gsea.1552686798457%5Cedb&pvalue=0.005&qvalue=0.1&overlap=0.5&similaritymetric=OVERLAP&combinedconstant=0.5&expressionfile=C%3A%5CUsers%5CHatim%5CDesktop%5CJunior%5Csem+2%5CKarras+Lab%5CGSEA%5CKarrasEXP.gct.txt HTTP/1.1[\r][\n]"
DEBUG - http-outgoing-4 >> "Host: localhost:1234[\r][\n]"
DEBUG - http-outgoing-4 >> "Connection: Keep-Alive[\r][\n]"
DEBUG - http-outgoing-4 >> "User-Agent: Apache-HttpClient/4.5.2 (Java/1.8.0_201)[\r][\n]"
DEBUG - http-outgoing-4 >> "Accept-Encoding: gzip,deflate[\r][\n]"
DEBUG - http-outgoing-4 >> "[\r][\n]"
DEBUG - http-outgoing-4 << "HTTP/1.1 500 Internal Server Error[\r][\n]"
DEBUG - http-outgoing-4 << "Date: Sun, 17 Mar 2019 06:55:04 GMT[\r][\n]"
DEBUG - http-outgoing-4 << "Content-Type: text/plain[\r][\n]"
DEBUG - http-outgoing-4 << "Access-Control-Allow-Origin: *[\r][\n]"
DEBUG - http-outgoing-4 << "Access-Control-Allow-Headers: origin, content-type, accept, authorization[\r][\n]"
DEBUG - http-outgoing-4 << "Access-Control-Allow-Credentials: true[\r][\n]"
DEBUG - http-outgoing-4 << "Access-Control-Allow-Methods: GET, POST, PUT, DELETE, OPTIONS, HEAD[\r][\n]"
DEBUG - http-outgoing-4 << "Content-Length: 20[\r][\n]"
DEBUG - http-outgoing-4 << "Server: Jetty(9.4.10.v20180503)[\r][\n]"
DEBUG - http-outgoing-4 << "[\r][\n]"
DEBUG - http-outgoing-4 << "Failed: empty String"
DEBUG - http-outgoing-4 << HTTP/1.1 500 Internal Server Error
DEBUG - http-outgoing-4 << Date: Sun, 17 Mar 2019 06:55:04 GMT
DEBUG - http-outgoing-4 << Content-Type: text/plain
DEBUG - http-outgoing-4 << Access-Control-Allow-Origin: *
DEBUG - http-outgoing-4 << Access-Control-Allow-Headers: origin, content-type, accept, authorization
DEBUG - http-outgoing-4 << Access-Control-Allow-Credentials: true
DEBUG - http-outgoing-4 << Access-Control-Allow-Methods: GET, POST, PUT, DELETE, OPTIONS, HEAD
DEBUG - http-outgoing-4 << Content-Length: 20
DEBUG - http-outgoing-4 << Server: Jetty(9.4.10.v20180503)
DEBUG - Connection can be kept alive indefinitely
INFO - status:Internal Server Error
INFO - Unable to create Enrichment Map:Internal Server ErrorResponseEntityProxy{[Content-Type: text/plain,Content-Length: 20,Chunked: false]}
DEBUG - Connection manager is shutting down
DEBUG - http-outgoing-4: Close connection
DEBUG - Connection manager shut down

Could you please explain the steps to reproduce this error? What data did you load, the settings and steps to connect to Cytoscape... and what do you expect from us. To me it seems like we should debug the GSEA software in order to solve this problem...

@ham123 Could you kindly review the answer below? I believe this should be definitely be upvoted and "accepted". The site member is knowledgeable in this specific problem.

This is a tough one, because the only relevant lines from your GSEA debug log were these two:
DEBUG - http-outgoing-4 << "Failed: empty String"
DEBUG - http-outgoing-4 << HTTP/1.1 500 Internal Server Error

My recommendations to you at the time would've been

try closing all the programs, restarting them, possibly rebooting and restarting them, and even analyzing a completely different set of data to see if you get the same results, —or—
update to the latest Java 8, reinstall Cytoscape 3.7.whatever and the EnrichmentMap plugin from scratch, and re-download GSEA, —or—
lastly, try reinstalling everything on a completely different machine/OS and see if you get different results (e.g., try on a Mac or Linux if you had problems with Windows originally, or vice-versa)

I know that's lame, but sometimes that's how you solve weird computer errors, when the error message the application itself gives you is non-specific and unhelpful.
If you had seen the error message
None of the gene sets have passed the filter. Try relaxing the gene set filter parameters.

in the GSEA application logs, then you would've had the problem described in EnrichmentMap's GitHub issue #358, and the solution would be to use less stringent p- and FDR Q-value cutoffs, as prescribed in the "Tips on Parameter Choice" section of the EnrichmentMap documentation.
My problem was that GSEA's default cutoffs (p-value < 0.005 and FDR Q-value < 0.1) were entirely too stringent for my results set, and using the "Very permissive" settings of p-value < 0.05 and FDR < 0.25 from the documentation solved it for me.
Hopefully this benefits someone who ends up here after googling for "Unable to create Enrichment Map:Internal Server ErrorResponseEntityProxy," which is just a stupendously unhelpful error message when "None of the gene sets have passed the filter" is actually what went wrong. Maybe go thumbs-up this comment in the EnrichmentMap issue tracker, if you think the maintainers need to work on that.


text:
Retrieving Metabolic networks from WIT Database
How can I get the dataset for all the 43 metabolic networks for this paper:
https://www.nature.com/articles/35036627
It looks like the WIT database does not exist anymore.

It seems like there are some files related to the databases at the ftp site: http://ftp.mcs.anl.gov/pub/WIT2/. The related files were last updated on 2002.
I found it via searching at google: site:*.mcs.anl.gov/ WIT.
Which reports also a page about the database. There seems to be a project to develop WIT3, because there is a mailing list for developers, maybe you can ask them more details if you have problems with the database.


text:
How do I tell featureCounts to ignore fragments that overlap with more than one feature?
I am using featureCounts from the subread package to count the number of fragments per gene in an RNA-seq experiment. I feed it a BAM file with paired-end reads aligned to a reference genome. In the help page of the executable, which I have installed using conda install subread, there is a section about overlaps, including for example the -O option which will assign the reads/fragments to all overlapping features. What is the default behaviour? Assigning to reads/fragments to the feature with largest overlap? not assigning multi-overlapping reads/fragments?
I would appreciate your help on this.
OS: macOS X
featureCounts version: 1.6.3

If a fragment overlaps multiple genomic features then it's unassignable, which is the appropriate way to handle ambiguous mappings when the resulting counts are then given to tools like DESeq2.

fair enough, so I will not use the -O option, thanks


text:
Splitting fasta file into smaller files based on header pattern
I have to split this fasta files into smaller files and write them into individual files my files
>lcl|CP000522.1_prot_ABO13860.1_1 [locus_tag=A1S_3471] [protein=hypothetical protein] [protein_id=ABO13860.1] [location=1..957] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13850.1_2 [locus_tag=A1S_3461] [protein=DNA replication protein] [protein_id=ABO13850.1] [location=950..1504] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13851.1_3 [locus_tag=A1S_3462] [protein=hypothetical protein] [protein_id=ABO13851.1] [location=complement(2523..3437)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13852.1_4 [locus_tag=A1S_3463] [protein=YPPCP.09C-like protein] [protein_id=ABO13852.1] [location=3538..4788] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13853.1_5 [locus_tag=A1S_3464] [protein=Cro-like protein] [protein_id=ABO13853.1] [location=5039..5629] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13854.1_6 [locus_tag=A1S_3465] [protein=hypothetical protein] [protein_id=ABO13854.1] [location=complement(6340..6906)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13855.1_7 [locus_tag=A1S_3466] [protein=Resolvase] [protein_id=ABO13855.1] [location=complement(7074..7685)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13856.1_8 [locus_tag=A1S_3467] [protein=hypothetical protein] [protein_id=ABO13856.1] [location=complement(8602..9732)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13857.1_9 [locus_tag=A1S_3468] [protein=putative lipoprotein] [protein_id=ABO13857.1] [location=complement(10072..10374)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13858.1_10 [locus_tag=A1S_3469] [protein=Diaminopimelate decarboxylase] [protein_id=ABO13858.1] [location=complement(10367..10723)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13859.1_11 [locus_tag=A1S_3470] [protein=regulatory protein LysR] [protein_id=ABO13859.1] [location=complement(12076..12444)] [gbkey=CDS]

The other pattern is
>lcl|CP000523.1_prot_ABO13861.1_1 [locus_tag=A1S_3472] [protein=DNA replication protein] [protein_id=ABO13861.1] [location=1..951] [gbkey=CDS]
>lcl|CP000523.1_prot_ABO13862.1_2 [locus_tag=A1S_3473] [protein=hypothetical protein] [protein_id=ABO13862.1] [location=3048..4262] [gbkey=CDS]
>lcl|CP000523.1_prot_ABO13863.1_3 [locus_tag=A1S_3474] [protein=hypothetical protein] [protein_id=ABO13863.1] [location=4357..5133] [gbkey=CDS]
>lcl|CP000523.1_prot_ABO13864.1_4 [locus_tag=A1S_3475] [protein=hypothetical protein] [protein_id=ABO13864.1] [location=6197..8608] [gbkey=CDS]
>lcl|CP000523.1_prot_ABO13865.1_5 [locus_tag=A1S_3476] [protein=secretory lipase] [protein_id=ABO13865.1] [location=8705..9403] [gbkey=CDS]

So now my idea is how do i parse and write them into individual files such that CP000522 output written to one file and CP000523 written to another file so forth and so on.
So far what i understand is i have to match the pattern after >lcl
so there are other patterns like "LN997847" in the file
Not sure how to proceed tried it in R but failed
it can be done with awk and sed which i tried but i can;t define something that parse all header like takes into account CP as well as LN .
Any suggestion or help would be highly appreciated .
here is my file

This is very easy with a bit of scripting one script will do it all. We do have R experts here (best solution)....Otherwise I'll post the code. How do you want each file named (important no matter which approach you take)?

Here's a simple awk approach:
awk '{if(/^>/){split($1,a,"[|.]")}print >> a[2]".fa"}' Protein_FASTA.txt

Or, more concisely, just:
awk '/^>/{split($1,a,"[|.]")}{print >> a[2]".fa"}' Protein_FASTA.txt

When run on the file linked to in your question, that results in the following files:
$ ls
AP014650.fa  CP003848.fa  CP007713.fa  CP012005.fa  CP015122.fa  CP017645.fa  CP018422.fa  CP020594.fa  CP023021.fa  CP024577.fa  CP026712.fa  CP027245.fa  CP030108.fa
CP000522.fa  CP003850.fa  CP007714.fa  CP012007.fa  CP015365.fa  CP017647.fa  CP018678.fa  CP020596.fa  CP023023.fa  CP024578.fa  CP026748.fa  CP027529.fa  CP030109.fa
CP000523.fa  CP003887.fa  CP008707.fa  CP012008.fa  CP015366.fa  CP017649.fa  CP018679.fa  CP021322.fa  CP023024.fa  CP025267.fa  CP026749.fa  CP027531.fa  CU459137.fa
CP000864.fa  CP003888.fa  CP008708.fa  CP012953.fa  CP015484.fa  CP017651.fa  CP019218.fa  CP021327.fa  CP023025.fa  CP026126.fa  CP027121.fa  CP027532.fa  CU459138.fa
CP000865.fa  CP003907.fa  CP008709.fa  CP012954.fa  CP015485.fa  CP017653.fa  CP020573.fa  CP021348.fa  CP023027.fa  CP026127.fa  CP027122.fa  CP027608.fa  CU459139.fa
CP001183.fa  CP003908.fa  CP008850.fa  CP012955.fa  CP015486.fa  CP017655.fa  CP020575.fa  CP021783.fa  CP023028.fa  CP026128.fa  CP027124.fa  CP027609.fa  CU459140.fa
CP001922.fa  CP003968.fa  CP008851.fa  CP012956.fa  CP016296.fa  CP017657.fa  CP020576.fa  CP021784.fa  CP023030.fa  CP026129.fa  CP027179.fa  CP027610.fa  JN377410.fa
CP001923.fa  CP004359.fa  CP010398.fa  CP013925.fa  CP016297.fa  CP018144.fa  CP020577.fa  CP021785.fa  CP023032.fa  CP026339.fa  CP027180.fa  CP029570.fa  LN865144.fa
CP001938.fa  CP006769.fa  CP010399.fa  CP014216.fa  CP016299.fa  CP018255.fa  CP020580.fa  CP021786.fa  CP023033.fa  CP026340.fa  CP027181.fa  CP029571.fa  LN997847.fa
CP002523.fa  CP007578.fa  CP010400.fa  CP014217.fa  CP016301.fa  CP018257.fa  CP020585.fa  CP021787.fa  CP023035.fa  CP026705.fa  CP027182.fa  CP029572.fa  LT594096.fa
CP002524.fa  CP007579.fa  CP010780.fa  CP014292.fa  CP016302.fa  CP018333.fa  CP020589.fa  CP022284.fa  CP024125.fa  CP026706.fa  CP027243.fa  CP029573.fa  Protein_FASTA.txt
CP003501.fa  CP007580.fa  CP010782.fa  CP014293.fa  CP017643.fa  CP018334.fa  CP020593.fa  CP022285.fa  CP024419.fa  CP026708.fa  CP027244.fa  CP030107.fa

Explanation

if(/^>/){split($1,a,"[|.]") : if this line starts with a >, split the first field on any occurrence of either | or . and save the results in the array a. Since your header lines all start with >lcl|, then the string you are looking for and a ., this means that the second value in the a array will be your target string.
print >> a[2]".fa" : print (append, >>) the current line to a file called "whatever the name of this sequence is" (a[2]) and .fa. This is run for every line in your input file. Note that if you run the same command again, you will need to first delete the files created the first time. If you don't, because I am using the >>, you will just append to the existing files.

I thought about perl -e  one-liner ... but I don't think a multiple file output is doable. The awk one-liner is cool.

Personally, I would have used a "_" along with "." and "|", but who am I. The small advantage of a script is that if its rerun it will clobber the old run output, without fear of appending.

@MichaelG. you can do it easily enough in perl, it's just slightly more cumbersome since you need to explicitly open and close each file: perl -ne 'if(/>.+?\|(.*?)\./){$name=$1; }open(my $fh, ">>","$name.fa"); print $fh "$_"; close($fh)' Protein_FASTA.txt. I didn't include the _ because the OP seemed to want the text before the version (the .N) to be the file name. I'd have used the _ myself too.

Thanks! I had no idea a file handle could be opened from "perl -e" ... I prefer awk code for this one.

wow so many solutions...let me run this..

thanks for the detailed explanation each line of the code..really appreciate it

You're welcome :). @MichaelG. note that I had a mistake in the first version of that perl one-liner. You need >>, not > (I've edited the comment) and it still has the same issue with existing files. The only way around that is to store everything in memory and only write at the end: perl -ne 'if(/>.+?\|(.*?)\./){$name=$1;}push @{$k{$name}},$_; END{for $name (keys(%k)){open(my $fh, ">","$name.fa"); print $fh @{$k{$name}}; close($fh)}}' Protein_FASTA.txt

I suggest replacing the if (…){…} block with the more awk-ish /^>/ {split($1, a, "[|.]")}, which simplified the actual print look to just {print >> a[2] ".fa"}.

@KonradRudolph that's how I would do it, yes. I used this version since I thought it would be easier for the OP to understand since it's more explicit and only 4 characters longer. But OK, I added the shorter version too.

@terdon I commented because I find the shorter version more readable (and I would argue that it is objectively so). Admittedly it does require knowing a bit about awk but that shouldn’t be an impediment.

@krushnachChandra ... its a culture. Make a code, condense a code and be more perlish, awkish or pythonic (as above).

I agree with @terdon flow control within one-line is minimised outside awk. This  awkish-ness so takes time to grasp.


text:
How to retrieve the best-scoring Trinity isoform from blastx results
I'm working on a somewhat unusual transcriptome focused on a killer X chromosome. This chromosome has some new genes, with expression levels quite different from their autosomal paralogs.
I've found that Trinity does an ok job of assembling the distinct genes, but when I try to do a gene-level DE analysis, it lumps the distinct genes together. So, I lose important information by doing a gene-level DE analysis and need to use a transcript-level approach.
I've done blastx with all of my transcripts using diamond, and want to extract the isoform from each gene with the best score (bitscore). I'm trying to write a script to parse this, but my scripting skills aren't up to the task. Does anyone know of (or have) a script that can do this?

Welcome to the site Ryan. Could you post that script so we can help you building these skills? Also it would be helpful if you explain how do you have organized those files/data

Hi Ryan! You mention that you lose important information by doing gene-level DE analysis. If you retain only a single isoform from each gene but still include isoforms from multiple paralogs, these will presumably continue to be lumped together, correct? On the other hand, if you retain only a single isoform from a single paralog, you've thrown away most of the interesting information. Perhaps you can clarify what you hope the result to be.


text:
How to interpret amino acid representation
What is the difference between C383S and C383A? For C383S, does it mean Cysteine(C) position 383 in the sequence and its next amino acid connect to it is Serine(S)?

This means that at position 383 the Cysteine(C) is mutated to a Serine(S).  Depending on the organism and gene usually the mutation notation is linked to a specific transcript since the position could vary depending on transcript splicing.

How about Pro136? Does it means Proline is deleted at position 136?


text:
Is there a safe catch-all adapter sequence for trimming?
I would like to trim/mark adapters using trimmomatic or picard MarkIlluminaAdapters from a series of Illumina Paired-End read fastqs. The fastq files may have been done using different kits or different generations/versions of illumina chips.
Without knowing ahead of time how the libraries were prepared (or which kit was used), Is there a safe "catch all" adapter sequence (or set of such sequences) that can be passed to a trimming tool and that would cover all Illumina kits? In other words, what are we expected to know about the fastq to ensure correct trimming?
edit: Picard defines a few built-in adapters which cover many standard kits, but it's not immediately obvious if you need to pass a common one (built-in sequences defined here).

You're best off just using fastp or Trim Galore!, both of which will determine the adapter sequence for you. Trim Galore! uses a built-in list of known sequences for this, whereas fastp uses read overlap.

For clarity, are you suggesting it's better to switch over completely to this recommended set of tools, or are you simply introducing the possibility of using fastp to determine the adapters (or extract the list from trim galore)?
One benefit of using MarkIlluminaAdapters is that it simply "marks" the sequences using XT:i: SAM trimming tag, rather than removing them altogether from the output bam/sam.

I see no benefit to using MarkIlluminaAdapters, replace it completely.

This article describes a "lossless" alignment pipeline which allows producing aligned bams which can be reversed back to the original fastqs. This is possible because of those trimming annotations, which annotate, rather than outright clip, the original sequence to indicate the portion of the read with adapters. I wish fastp supported these SAM trimming annotations, it would be a strict improvement.

If you want that just use an aligner with soft clipping. Note also that few aligners support SAM input.

Lol, I think I'll start calling mappers "read slingers"  To be serious though, it's possible for a conscientious bioinformatician to store reads in SAM (et al.) format at every stage (raw unaligned, QC'd, aligned, aligned post-processed, etc.) in a way that preserves the original raw sequence unambiguously. In practice though, it will be problematic since many tools are not designed with this paradigm in mind. It would take years to get everyone on the same page, if indeed it's even possible to do so.

Lol, stupid autocorrect on my phone.

re: aligners not supporting SAM. small note. The workflow I mentioned above will temporarily transform the unaligned SAM (with XT:i adapter trimming annotations) into fastq that has adapters hard-clipped. The mapper is run on that clipped fastq (like you say, because most of them only support fastq). And then this tool, MergeBamAlignment, combines the mapping information from the aligned BAM into the full, soft-clipped, reads in the unaligned SAM. It might be straighforward to compute the XT:i annotations by comparing the fastp input/output , but I haven't fiddled with it. I'll consider.

Alternatively, I really like using bbduk which is part of the BBMap suite.
I've processed every nascent sequencing dataset that has been published, and found a lot of quirky errors with older datasets using TrimGalore. bbduk is a little more fine-tuneable relative to cutadapt/trimmomatic/trimGalore (built on top of cutadapt)/fastp and the run time and memory usage is significantly less.
When you clone/install the repository, there is a directory of adapter fa files that you can choose from for trimming in the resources/ subdirectory which you can specify using the "ref" argument (i.e. ref=adapters.fa). The "adapters.fa" file is an excellent catchall for adapters, but you will likely over-trim a little. Alternatively, if you know what adapters were used in your dataset, you can choose one of the other pre-configured fa files. There are also a number of stats/marking/reporting/trimmed read save options which you can see here:
https://github.com/BioInfoTools/BBMap/blob/master/sh/bbduk.sh


text:
How to interpret UniProt allele patterns?
In creating a parser for the UniProt flat file uniprot_sprot.dat (directory) in the manual for the comments POLYMORPHISM there are listed patterns for alleles, e.g. for 1A11_HUMAN is

The following alleles of A-11 are known:
  A*11:01 (A-11E), A*11:02 (A-11K), A*11:03, A*11:04, A*11:05 and A*11:07.

I tried searching for a description of the patterns using the words allele and/or polymorphism and found nothing.

How are the allele patterns to be interpreted?
Is there a search term to be used when looking for the patterns?
Also, is there an official repository of the alleles?

The answer does not have to explain the specifics of the example given.

Example of explanation of other allele pattern from here

This explanation of the patterns doesn't explain the part of the pattern in UniProt like (A-11E).

Of interest: Nomenclature for Factors of the HLA System

These are cytotoxic T-cell HLA alleles. HLA genotyping is very common and easy to do, so Genbank is the repository.
A*11:01 has very high frequency in Aborigine populations here . You can explore the population genetics and the past, present and current population genetics distributions per population of HLA at http://www.allelefrequencies.net/

Thanks. I gave it an up-vote, but I don't see how this answers my main question of How are the allele patterns to be interpreted? Also my background is in computer programming, so please don't take for granted I have knowledge on the bioinformatic information. :)

Thanks, I would look at whether these alleles cluster geographically. Secondly, the diseases which induce high HLA*11:  alleles

My goal is to try and parse the allele patterns further than what they currently are (strings of characters) so that the parser can validate that they are correct. Without the pattern the only way to verify them would be to check each one against an database. Even then the allele might be correct but not in the database. With a pattern there is a second means of validating the allele.

HLA alleles should have been extensively characterised. I have a colleague that types this stuff alot, but I tend to avoid it ( complicated ethics applications)


text:
Preparing PDB file for modeling with Swiss-Model
I recently starting studying bioinformatics and have a question. I need to build a model of protein using Swiss-Model. First, I need to download a protein record from PDB (for example, 2HZI), then I must open PDB file using simple text editor. I know this usually requires biopython, but i don't use it yet. Now I'm trying to study how to build a model.
Swiss-model needs only one chain to build model, and I remove from file all extra information like AUTHOR,REMARK,HETATM etc. But Swiss-model still gives an error: "Sequence String contains illegal characters. Only letters and gaps allowed."

Welcome to the site @Badscience. Perhaps you haven't noticed but the site language is English. If you could translate your question into English we might be able to help. Please [edit] the question in English. About which program do you have a problem? I see Biopython and PDB mentioned but you might be trying to parse a file or using it to connect with another program.

Thank you very much, @llrs. I'm sorry if i do some mistakes in this text.
I am just starting to study bioinformatics and i have a problem.  Today i download protein from PDB and try to build a model, using "Swiss-model" site, but this is requires work with PDB file, i don't use biopython yet, because i only try to understand how to build a model. I should delete  excess information from PDB file, information like MODEL.ENDMDL,REMARK etc. But when i try to prepare file to modeling, i still face error: "Sequence String contains illegal characters. Only letters and gaps allowed".

please [edit] your question with the details of the comment and more: Which "Swiss-model" site? (add a url at least). What are you trying to build? A model of a protein, or a complex or something else? Could you add the text of the file you are uploading to that site? Also don't worry much about the English, we would help you improving it if we can.

As explained in the website the data introduced must be:

Format must be FASTA, Clustal, plain string, or a valid UniProtKB AC

So you need to paste just the protein sequence. Without seeing the data you paste I would guess that you don't use a valid plain string.
As Jonathan Moore explains below a FASTA format is:
>SEQUENCE_1
MTEITAAMVKELRESTGAGMMDCKNALSETNGDFDKAVQLLREKGLGKAAKKADRLAAEG
LVSVKVSDDFTIAAMRPSYLSYEDLDMTFVENEYKALVAELEKENEERRRLKDPNKPEHK
IPQFASRKQLSDAILKEAEEKIKEELKAQGKPEKIWDNIIPGKMNSFIADNSQLDSKLTL
MGQFYVMDDKKTVEQVIAEKEKEFGGKIKIVEFICFEVGEGLEKKTEDFAAEVAAQL

(breaking the sequence in several lines is optional but usually done fore easy reading)

I agree with this comment.  Take a look at the FASTA format and present your protein in a text file in that format, like this:
{>SEQUENCE_1
MTEITAAMVKELRESTGAGMMDCKNALSETNGDFDKAVQLLREKGLGKAAKKADRLAAEG
LVSVKVSDDFTIAAMRPSYLSYEDLDMTFVENEYKALVAELEKENEERRRLKDPNKPEHK
IPQFASRKQLSDAILKEAEEKIKEELKAQGKPEKIWDNIIPGKMNSFIADNSQLDSKLTL
MGQFYVMDDKKTVEQVIAEKEKEFGGKIKIVEFICFEVGEGLEKKTEDFAAEVAAQL}

What they should be doing is overlaying the Swiss-Prot prediction against the real sequence data .... thats when things get harder


text:
What is the purpose of the "fixedSearchFolder" and "fixedCombinedFolder" in MaxQuant?
What is the purpose of changing the fixedSearchFolder and fixedCombinedFolder in MaxQuant proteomics software? Is there any data that can be reused in other runs, so would it be good to use the same search folder for multiple runs? Or is it simply to direct temporary data and results data into different locations and each run should use a different folder?


text:
Matching these matrices in R
I have two matrices; I want to convert the row names of first matrix to gene symbol from matched ensemble=gene symbol from second matrix
> head(mat1[,1:2])
                TCGA-L5-A4OG-11A-12R-A260-31
ENSG00000000003                         1818
ENSG00000000005                            0
ENSG00000000419                         1436
ENSG00000000457                         1175
ENSG00000000460                          242
ENSG00000000938                          536
                TCGA-IC-A6RE-11A-12R-A336-31
ENSG00000000003                         4596
ENSG00000000005                            3
ENSG00000000419                          751
ENSG00000000457                          840
ENSG00000000460                          205
ENSG00000000938                          253
>

> dim(mat1)
[1] 56925    11
>

> head(mat2)
                ensembl              symbol
ENSG00000274572 ENSG00000274572       ZYXP1
ENSG00000159840 ENSG00000159840         ZYX
ENSG00000162378 ENSG00000162378      ZYG11B
ENSG00000232242 ENSG00000232242    ZYG11AP1
ENSG00000203995 ENSG00000203995      ZYG11A
ENSG00000070476 ENSG00000070476        ZXDC
>

> dim(mat2)
[1] 36848     2
>

How I can do that in R?
Thank you

I'm voting to close this question as off-topic because it is not about bioinformatics

Since they share the ensembl_ID column, you can merge them, then assign the symbol column to the rownames, then delete the symbol column.
something like:
merged <- merge(mat1, mat2, by = 0)
rownames(merge) <- merged[,'symbol']
merged[,symbol] <- NULL

https://stat.ethz.ch/R-manual/R-devel/library/base/html/merge.html


text:
Pipeline for extracting gene from multiple genomes for use in HyPhy selection analyses?
I have been trying to obtain some preliminary data from HyPhy selection analyses to inform a larger project. I have obtained a number of assembled mammalian genomes from NCBI with the initial goal of extracting from each of them their corresponding sequence for a specific gene. I wanted to use HyPhy to determine if there are signatures of positive selection in any of the mammalian lineages for that specific gene. However, I have been struggling to make headway with this for some time and am realizing that I may be out of my depth.
I was hoping that people could suggest, or at least point me in the direction of, a pipeline/guide/best practices for how one would go about achieving what I am trying to do. I would like to start from scratch and make sure I am performing my analyses in a proper way, using advice from other users.
As I mentioned, I currently have a collection of assembled genomes from NCBI which I would like to pull specific genes from. My understanding is that I will then need to extract the exons (I imagine there is a way to do this without having to extract the entire gene first), align the sequences from each species while preserving reading frames, and then use that alignment as input in HyPhy's codon models. It is very possible that I am missing out some fundamental steps in the process.
If anybody could offer direction/advice, I would be very appreciative.

Very interesting, personally I'd use PAML. @NatWH is good advice. Keep us posted.

That is interesting. I believe that others in my lab have more experience with HyPhy, hence why that is the goal at the moment. I will look into PAML as well. Thank you!

Okay well, because you're using HYPHY, you're going to need a gene tree of every gene you intend to analyze. This requires you to understand the homology relationships between the genomes at hand. Depending on your project, you might already have a list of homologs you're trying to analyze, or you can start from scratch.
If starting from scratch, you will need the proteomes of the species you're interested in (that is, the translation of the coding sequences of each gene, or to put it another way, the concatenated exons). Genome releases should typically have a proteome available for download, along with its corresponding coding sequence.
Next, you need some way to get the homolog groups for the genes in your genomes. You can use a pipeline like OrthoFinder, or iterative tree building and trimming, to infer homolog clusters (OrthoFinder calls them orthogroups but they're the same thing - homologous genes for the species of interest including both single copy orthologs and paralogs).

For this step, if you already have a set of genes and you know the copy numbers, or know that they're single copy orthologs, then it is more easy - you can just pull the sequence from each proteome and go from there.

Once you have your homologs, you need to align the amino acids using a tool like MAFFT, and then align the coding sequences to the amino acid alignment to attain a coding sequence alignment. For this you can use a tool like pal2nal, or check out pxaa2cdn from the phyx package written by some of my collaborators.

If the sequences aren't too divergent, and intron-exon structure is conserved, you could possibly also include alignments of the introns to provide more phylogenetic information. If you do use this approach however, you'd need to align the exons and introns separately. Also, the introns would not be included in the selection analysis.

Next, I would typically reinfer the gene tree of each gene using nucleotide models, which depending on the time, level of selective constraint, etc. can give more signal for the phylogenetic relationships. For this you can use a tree inference program like RAxML. With this tree and your coding sequence alignment, you'll be ready to run analyses in HYPHY (aBS-REL, MEME, FUBAR, etc., depending on your requirements).

Note that if you want, you can also do this in a Bayesian framework to integrate over uncertainty in the model parameters, tree topology, etc. Several Bayesian phylogenetic inference tools (e.g. MrBayes, PhyloBayes), allow you to infer putatively positively selected sites at the same time as tree inference.

Thank you, NatWH! Your answer seems perfect - comprehensive and easy to follow! I will try to follow this pipeline.

@PollardMD no problem, feel free to comment here if you run into any problems.


text:
Missing lower IQR in boxplot
Im not able to figure it out why the lower quartile of the two data points missing.
Data
Code
require(reshape2)
require(ggplot2)
require(RColorBrewer)
#require(dplyr)
require(tidyr)
require(tibble)
df <- read.csv("NC_II_box_plot_input.txt",sep = ",")

df_melt = melt(df,id.vars = NULL)

df_melt$Group <- gsub("HSC",'',df_melt$variable)
head(df_melt)
df_melt$Group <- gsub("CMP",'',df_melt$Group)
df_melt$Group <- gsub("GMP",'', df_melt$Group)
df_melt$Group <- gsub("Mono",'', df_melt$Group)
df_melt$Group <- factor(df_melt$Group,levels = c("m","nc","r","rr","sn","t"))

ggplot(df_melt, aes(variable,value)) +
  stat_boxplot(geom="errorbar", width=.5)+
  geom_boxplot(aes(fill=Group), position = position_dodge2(preserve = "total"))+
  scale_fill_manual(values=c("grey","dodgerblue4","red","darkgreen","brown","maroon1")) + #for boxplot
  # scale_fill_manual(values=c("grey","dodgerblue4","red","darkgreen","brown")) + #for boxplot
  geom_jitter(position=position_jitter(width=0.3), size=0.01, colour="black") +

  #scale_fill_manual(values = c("black","brown","maroon1","darkgreen"))+
  #stat_summary(fun.y=mean, colour="red", geom="line", aes(group = 1))+
  #geom_jitter(position = position_jitter(0.0001),size=.00001)+
  #theme_bw(base_size=15)+
  theme_bw()+

  theme(axis.text.x=element_text(angle = 45, size=25, face="bold", hjust = 1),
        axis.text.y=element_text(angle=0, size=35, face="bold", vjust=0.5),
        plot.title = element_text(size=40, face="bold"),
        legend.title=element_blank(),
        legend.key.size=unit(1, "cm"),      #Sets overall area/size of the legend
        legend.text=element_text(size=30)) +

  labs(x="Cell_Type", y=expression(paste("Expression"), title="Expression")) +
  guides(colour = guide_legend(override.aes = list(size=5)))

#

to find the quartile i used this as well
A <- boxplot(df)
mytable <- A$stats
colnames(mytable)<-A$names
rownames(mytable)<-c('min','lower quartile','median','upper quartile','max')
mytable

Any help or suggestion would be highly appreciated

Could you clarify what is the quartile you can't find? (making the example more minimal) Do you mean that you want it to be plotted and it isn't or something else?

if you can look then the rrCMP lower quartile  is missing is it because the data is as such or due to some code error that im not sure..about ..

That is not the quantile but the 1.5*inter-quartile-range. Could you just provide the data for these category? Besides this question is purely about plotting/statistics not about bioinformatics (an expert on ggplot2 would be in much better position to answer if it is a bug on its side)

sorry yes the IQR my bad .Yes I have given the data file in the start of my question...

Well, to my taste this is too much data, your problem doesn't need all the data you provided. It will help to focus the question...


text:
Attractor Landscape Analysis
I have come across a modeling toolbox, ATLANTIS, which is able to determine cell fates in silico based on the input models provided. This MATLAB-based toolbox is built on a method called "Attractor Landscape Analysis." I tried to find a paper dedicated to explaining this method but I did not have any luck so far. Can any of you folks please explain it to me, especially the term "attractor"?

From Wikipedia:
In the mathematical field of dynamical systems, an attractor is a set of numerical values toward which a system tends to evolve, for a wide variety of starting conditions of the system.
In biology you can think of an attractor as a state towards the system is moving to. A good example is the process of differentiation where stem cells compromise to more specialized cells, that could be thought as the attractors of the process.

The problem I have is attractor theory was derived from chaos theory. If you didn't know that you should. You need to understand the 3 components of chaos theory. BTW, the word "evolve" is a total misuse of the term 'evolution' and is not even used in chaos theory.  If you are not prepared to understand chaos theory you will not understand attractor theory.

@MichaelG. can you please recommed a good resource to understand chaos theory in the context of biological sciences?

Any basic text book will summarise chaos theory, even very basic popular science (Dummies Guide to Chaos, type books). My information is rooted in Darwinian selection vs. chaos theory, and unless you are an evolutionist you wouldn't understand what the debate was about, who won the debate and why.

Oh sorry and ignore anything on fractals, it's not relevant.

@MichaelG. this video explains Chaos and Attractor Theory beautifully. https://www.youtube.com/watch?v=cwaRtew-CIU . However, I might have confused myself with the definitions of "internal noise" and "external perturbations." In the context of biology, can you please explain the difference between the two?

Thanks. Let me summarise, because we're a bit off-topic, this is more like a Stats Exchange question. Attractor theory risks parameter sets that will produce weird results. The weird results are chaotic artefacts. The explanation as to why is more complicated.


text:
mapping heteryzygous kmers on a genome
I have a set (a couple of millions) of kmer pairs of known length (usually 21) that I know that are heterozygous in the middle nucleotide. For instance:
AAAAAAAAAANAAAAAAAAAA

Where N has always two states. I would like to map these nucleotides on a fragmented genome assembly that is presumably haploid. I bet that this problem is resolved, but I have hard times to make a practical solution, preferably within python environment. I considered following

Using an existing mapper

I tried to use python binding of minimap2, thinking that I just need to tell it that I search for exact matches only, but I did not manage to set it up. I set up smaller minimisers but it always returned no hits for kmers of length 21 (I tried to tweak parameters k, w and min_cnt, but to be honest I don't have a good understanding of how minimizers work). Is there a mapper that would be easily set up to do these searches?

Hash table

The second thing I tried was to build a hash table of all k - 1 / 2 mers of the genome to map the left and right sides of the kmers. Where sequences were mapped to lists of corresponding genomic positions, but this have blown my memory very soon (~50m of kmers were indexed using ~8GB of memory; i.e. this approach does not seems to be viable on a genomic scale).
I am also suspecting that there will be an elegant solution using a suffix array, but I have not cracked that one yet.

I don't know if this is the best answer, but one quick-n-dirty approach I commonly use involves tweaking bwa mem to report only perfect matches. So in your case, the idea is to make two copies of each 21-mer, replace the N with the appropriate allele, and then map the resulting sequences to your bwa-indexed assembly like so.
bwa mem -k 21 -T 21 -a -c 5000 assembly.fasta your-21mers.fasta

The -k 21 says to use a minimum seed length of 21, which forces exact matches. The -T 21 requires a minimum score of 21, which also enforces an exact match. The -a parameter reports all matches, since a "best" match doesn't make sense in this situation. The -c parameter limits how many matches are reported, which may need to be adjusted depending on how repetitive the 21-mer is.
UPDATE If you don't want to introduce a dependency like bwapy, it's fairly straightforward to implement this in pure Python.
with tempfile.TemporaryFile() as kmerfile:
    for kmer in kmers:
        print('>{k}\n{k}'.format(k=kmer), file=kmerfile)
    kmerfile.flush()             # Often works without these two lines
    os.fsync(kmerfile.fileno())  # But best to be safe
    bwacmd = ['bwa', 'mem', '-k', ksize, '-t', ksize, '-a', '-c', '5000', refrfile, kmerfile.name]
    bwaproc = subprocess.Popen(bwacmd, stdout=subprocess.PIPE, universal_newlines=True)
    for line in bwaproc.stdout:
        # process SAM records line by line

Also, one can call bwa mem via a C API, but the mechanics of this are left as an exercise for the reader.

Indeed, with suffix arrays, you can map a kmer in ~0.0007s. The way to go is to build a suffix array (for instance using PySAIS package) and then do a convoluted binary search of the four possible kmers. A full implementation of such mapper can be found here.
To build suffix array
from Bio import SeqIO
from Bio.Seq import Seq
from PySAIS import sais
from bisect import bisect_right

# load the genome
kmer_genome_file = 'my_genome.fasta'
ffile = SeqIO.parse(kmer_genome_file, "fasta")
sequences = []
scf_names = []
for seq_record in ffile:
    scf_names.append(seq_record.name)
    sequences.append(str(seq_record.seq).upper())

ffile.close()

# paste all scaffolds into a single string and build a suffix array
genome = '$'.join(sequences)
sa = sais(genome)

To transform genome coordinate to scaffold/position coordinates we define a function evaluated2assembly_position
def evaluated2assembly_position(eval_pos):
    scf_index = bisect_right(scf_sizes, eval_pos)
    eval_pos - scf_sizes[scf_index - 1]
    if scf_index > 0:
        return((scf_names[scf_index], eval_pos - scf_sizes[scf_index - 1]))
    else:
        return((scf_names[scf_index], eval_pos))

and finally the function for a convoluted binary search (assuming all kmers being 21 nt long and having the unknown nucleotide in the middle, on 11th nt). Note that you need to map to both strands of DNA, which is equivalent of mapping kmer and it's reverse complementary sequence
def searchKmer(self, kmer):
    reverse_complementary_kmer = str(Seq(kmer).reverse_complement())
    alignments = []
    for strand in ['+', '-']:
        if strand == '-':
            kmer = reverse_complementary_kmer
        L_kmer = kmer[0:10]
        l = 0
        r = len(sa) - 1
        while l <= r:
            m = (l + r) // 2
            eval_pos = sa[m]
            L_genome_kmer = genome[eval_pos:(eval_pos + 10)]
            if L_kmer == L_genome_kmer:
                r_L_converged = r
                for N in ['A', 'C', 'G', 'T']:
                    kmer = kmer[0:10] + N + kmer[11:21]
                    while l <= r:
                        m = (l + r) // 2
                        eval_pos = sa[m]
                        genome_kmer = genome[eval_pos:(eval_pos + 21)]
                        if kmer == genome_kmer:
                            m_hit = m
                            # following ms
                            while kmer == genome_kmer:
                                scf, pos = evaluated2assembly_position(eval_pos)
                                alignments.append([scf, pos, strand, N])
                                m += 1
                                eval_pos = sa[m]
                                genome_kmer = genome[eval_pos:(eval_pos + 21)]
                            # previous ms
                            m = m_hit - 1
                            eval_pos = sa[m]
                            genome_kmer = genome[eval_pos:(eval_pos + 21)]
                            while kmer == genome_kmer:
                                scf, pos = evaluated2assembly_position(eval_pos)
                                alignments.append([scf, pos, strand, N])
                                m -= 1
                                eval_pos = sa[m]
                                genome_kmer = genome[eval_pos:(eval_pos + 21)]
                            break
                        elif genome_kmer < kmer:
                            l = m + 1
                        else:
                            r = m - 1
                    r = r_L_converged # recover the r after converging on the left sub-kmer
                break
            elif L_genome_kmer < L_kmer:
                l = m + 1
            else:
                r = m - 1
    return(alignments)

Taking one of the fist kmers in the scaffold one to test it
searchKmer('ACGGCCACCCNGTGTCGTTGT')
Out[23]: [['scaf000001', 28, 'G']]

And a softmasked kmer that we can expect more than once
kmer = 'ccgaaccacgNtggaggcctg'.upper()

searchKmer(kmer)
Out[31]:
[['CCGAACCACGATGGAGGCCTG', '1_Tps_b3v07_scaf000306', 5196, 'A'],
 ['CCGAACCACGCTGGAGGCCTG', '1_Tps_b3v07_scaf008224', 15781, 'C'],
 ['CCGAACCACGGTGGAGGCCTG', '1_Tps_b3v07_scaf001152', 71261, 'G'],
 ['CCGAACCACGTTGGAGGCCTG', '1_Tps_b3v07_scaf002780', 1721, 'T']]

So to map all your kmers (assuming a file with a kmer per line)
kmer_file_name = 'my_kmer_file.txt'
with open(kmer_file_name, 'r') as kmer_file:
    kmers = [kmer.rstrip() for kmer in kmer_file]

mapping_list = [searchKmer(kmer) for kmer in kmers]

Stuff can be probably done in a more elegant way, but this is a working solution that scales to any genome size.


text:
How snippy makes MSA-like aligned fasta output from multiple samples?
From the log file it seems snippy doesn't do assembly. It only does mapping:
fastq --> SAM --> BAM --> VCF --> consensus_seq/snps

But if multiple samples are given (with snippy-multi option) it outputs an aligned fasta format. My understanding so far is it doesn't go through any MSA alignment algorithm, but use the pairwise alignment information from SAM files.
How is the aligned fasta file generated? The log file doesn't contain information about how snippy calculates it.  Is there any other tool that I can use to get similar aligned fasta file from SAM/VCF file without going through any actual MSA algorithm?

I tried to make the question bit more clear, but I am not sure if I got right what you were asking for (maybe that's why I had a need to try to improve it), please, don't hesitate to [edit] the question if I the changes I made were inappropriate.


text:
log2fold change as heatmap
From a paper:

(D) Expression analysis of multiple lineage-specific differentiation markers in WT and PUS7-KO EBs (14 days). Heatmap shows log2 fold change (FC) PUS7-KO to WT for each individual gene (rows) in three independent experiments (columns).

They have analysed the data in EdgeR but I was wondering how did they plot fold change when they have two condition but they have three independent experiment, I haven't used EdgeR so far, use DESeq2 where I group my biological replicate. For example if I have two condition with three biological replicate then I would group them and analyse, in case of DESeq2 I would get final result with one fold change report.
How did they do in this case reporting three fold change ?
Data sets where I see three wild type and three knockout.

If I read quickly thru the paper, I see you refer to Fig 5D. In the methods they describe that they use Fluidigm (i.e. RT-qPCR) for gene expression, not RNA-seq. The data link you give is for small-RNA-seq.

yes fig5D ."In the methods they describe that they use Fluidigm (i.e. RT-qPCR) for gene expression, not RNA-seq." okay i think i have to read it ,i missed it .


text:
Old versions of the reference genome and dbSNP
For benchmark purposes I’m trying to find an old version of the human reference genome (pre-GRCh37) and a matching version of dbSNP.
Unfortunately it appears as though older versions of dbSNP aren’t archived, or at least these archives are not publicly accessible.
In particular, the NCBI FTP archive is fairly sparse. The only pre-GRCh37 version it contains is ostensibly dbSNP build 124 from January 2005, for hg17 (NCBI build 35). Unfortunately the corresponding FTP folder is actually empty. dbSNP build 138 has been “removed to recover storage”, whatever that means, and I’m unclear which version of the human reference genome it refers to (its official page does not contain this rather crucial bit of information).
Is there another authoritative source of this data?

Perhaps ENSEMBL will have some of what you need? At least the older reference genomes appear to be there. https://www.ensembl.org/info/data/ftp/index.html

@Scot Ensembl has variation data but it does not, as far as I know, host the dbSNP data (i.e. a dataset exactly corresponding to a given dbSNP build). Although I might be able to filter the corresponding entries from their variant database, and convert it to VCF.

Actually this won’t work: The NCBI36 variant data on Ensembl doesn’t contain dbSNP information directly, and it’s in a format that’s not trivial to convert to VCF.

The UCSC Genome Browser downloads page has old human genome assembly versions way back to the first human genome ever published (by UCSC in 2000).
http://hgdownload.soe.ucsc.edu/downloads.html#human
We also have a variety of archived dbSNP data up to dbSNP141 on hg38, dbSNP 131 on hg37, dbSNP126 on hg18:
http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/snp141.txt.gz
http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/snp138.txt.gz
http://hgdownload.soe.ucsc.edu/goldenPath/hg18/database/snp126.txt.gz
They're in a custom table format, which can be turned into VCF manually.
You can also get all of this data from Table Browser and add filters, etc.:
http://genome.ucsc.edu/cgi-bin/hgTables


text:
How can I convert Affymetrix microarray data to bedgraph format?
I've previously asked a different question about the ChIP-chip microarray data that I'm working on. I'm new to this so I'm sorry for asking here again.
I have .CEL files and I'm using the Affy package for R to read them. After I load all the files in, I then normalise the data and  generate a .txt file containing the expression values of each file (ranging from 11 to 12). Column 1 contains the gene names and the rest of the  columns contain the expression values for each of the files.
I was wondering if there is a way to convert these data to bedgraph or BED format? What I'd like to do is a metagene analysis but I don't know how that would be possible with expression value data.
I looked online and there are some kinda relevant suggestions/solutions, but they are for Windows machines and I'm working on MacOS.

Do the generated txt files contain the physical coordinates of the probes? If so, converting them to BED should be a simple matter of rearranging the columns (or they might already be valid BED files).

Hi Konrad, thank you for your response.
The generated txt file does not contain the coordinates of the probes, however, I have a gff file that contains all this information and I guess I could look up the coordinates for each gene of the array. However, I don't know how the expression values could be converted?

Why makes the solutions for Windows specific for Windows? Did you try to overcome these using equivalent programs or formats?


text:
samtools view: writing to standard output failed: Broken pipe
I have the following alignment script which uses BWA:
#!/bin/bash
#usage: sh align_pbs.sh /scratch/waterhouse_team/benth/dbg2olc-40x/good/ output.racon2.fasta bam-soap

mkdir $3

for r1 in $(find $1 -name "*R1*.gz");
do
  base=${r1%_R1*}
  output=$(basename $(echo $r1 | sed 's/R1//g'))
  r2=$(echo $r1 | sed 's/R1/R2/g')

  #Get read group infomration:
  #Source: https://www.biostars.org/p/280837/#310132
  header=$(zcat $r1 | head -n 1)
  id=$(echo $header | head -n 1 | cut -f 1-4 -d":" | sed 's/@//' | sed 's/:/_/g')
  sm=$(echo $header | head -n 1 | grep -Eo "[ATGCN]+$")
  #echo "Read Group @RG\tID:$id\tSM:$id"_"$sm\tLB:$id"_"$sm\tPL:ILLUMINA"
  echo "@RG\tID:${output}\tSM:${output}\tLB:${output}\tPL:ILLUMINA"

  #cat <<EOF
  qsub <<EOF
#!/bin/bash -l

#PBS -N $output
#PBS -l walltime=48:00:00
#PBS -j oe
#PBS -l mem=70G
#PBS -l ncpus=12
#PBS -M m.lorenc@qut.edu.au
##PBS -m bea

cd \$PBS_O_WORKDIR

conda activate bwa-hic
bwa mem \
$2 $r1 $r2 \
-t 12 \
-R "$(echo "@RG\tID:${output}\tSM:${output}\tLB:${output}\tPL:ILLUMINA")" | samblaster -r | samtools view -@ 12 -h -f 0x2 -F 2316 - | samtools fixmate - - | sed 's/^/LP1-/' | samtools sort -@ 12 - -o ${3}/${output}.sorted.dedup.bam

EOF

done

Unfortunatelye, this pipeline cuased the following erros:
samblaster: Version 0.1.24
samblaster: Inputting from stdin
samblaster: Outputting to stdout
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 1200000 sequences (120000000 bp)...
[M::process] read 1200000 sequences (120000000 bp)...
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (2, 316829, 6, 4)
[M::mem_pestat] skip orientation FF as there are not enough pairs
[M::mem_pestat] analyzing insert size distribution for orientation FR...
[M::mem_pestat] (25, 50, 75) percentile: (432, 474, 509)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (278, 663)
[M::mem_pestat] mean and std.dev: (471.59, 53.66)
[M::mem_pestat] low and high boundaries for proper pairs: (201, 740)
[M::mem_pestat] skip orientation RF as there are not enough pairs
[M::mem_pestat] skip orientation RR as there are not enough pairs
[M::mem_process_seqs] Processed 1200000 reads in 615.287 CPU sec, 58.490 real sec
samblaster: Loaded 19 header sequence entries.
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
sed: couldn't write 567 items to stdout: Broken pipe
[E::bgzf_flush] File write failed (wrong size)
samtools fixmate: Couldn't write to output file: Broken pipe
[E::bgzf_close] File write failed
[bam_mating] error while closing output file
samtools view: writing to standard output failed: Broken pipe
samtools view: error closing standard output: -1
samblaster: Unable to write to output file.

Thank you in advance,

[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting

SAMtools sort has been unable to parse its input, which it thought was SAM (mostly because it couldn't be recognised as another format e.g. BAM).
This is because sed 's/^/LP1-/' is putting LP1- at the front of every line. If the output of samtools fixmate is SAM, then this LP1 is garbling the SAM header lines. If the output of samtools fixmate is BAM (as it probably is), this LP1 is completely garbling the data.
What were you trying to achieve with this sed command?

It is attended that `sed 's/^/LP1-/' is putting LP1- at the front of every line because I would like to use SGSautoSNP.py to detect different cultivars. How is it possible to fix the above pipeline?


text:
Voom transformation of RNA seq raw counts data
I have RNA seq data of raw gene counts that I want to transform for linear modelling. I am trying to voom transform, to do a weighted analysis.
Data frame: prac_count_10
gene         S1     S2    S3   S4    S5
ENG0000456   0      1     10   145   24
ENG0000458   7      2     9    0     0
ENG0000657   76     12    56   10    2
ENG0000689   0      0     0    3     5

Code:
prac_prac_10 <- voom(prac_count_10[,2:5], design=NULL, lib.size=NULL,
  normalize.method = "none", span = 0.5

However this doesn't give an output with the gene sample names and logCPM values. Would like an output with gene id's and logCPM values for linear modelling.
Also tried:
cpm <- cpm(prac_count_10)
lcpm <- cpm(prac_count_10, log=TRUE)

````````````

Try to use DElist() function before you transform, and also make rownames first.
counts <- prac_count_10[,-1]

rownames(counts) <- prac_count_10[,1]

DGE1 <- DGEList(counts)

And then continue with your voom...

Thank you. what package is DElist in?

That would be edgeR.

If I do:   prac_count <- data.frame(prac_count, row.names=1), counts<- prac_count[,1], rownames(counts) <- prac_count[,1], DGE1 <-DGEList(counts)  If doesn't work... have I missed a step? Thank you so much

That's not the code I gave in my answer, is it?

Don't mention it! Good luck with it.


text:
Are phylogenetic tree construction algorithms any different than general clustering algorithms?
Are phylogenetic tree construction algorithms any different from clustering algorhithms? I suspect the answer is no.
Of course phylogenetic tree construction uses biological knowledge, e.g special distance metrics, but does it brings anything new to clustering algorithm (e.g hierarchical, neighbor-joining etc).

Yes, they are very different. If no one else responds I'll provide a full theoretical perspective about what is going on.

@Michael G. Eagerly waiting for your answer. :)

Most phylogenetics theory is based on the probability of reversion mutations and this drastically distinguishes it from cluster algorithms based on a simple distance matrix

Well, distance-based tree building algorithms are based on a simple distance matrix.

The goal of a phylogeny is to estimate the "expected" number of mutations between all taxa in the analysis and their hypothetical common ancestors. A cluster-analysis will only identify the "observed" mutations and "expected" and "observed" mutations can be majorly different due to the major artefact of reversion mutation. This is particularly true of nucleotide phylogenies.
The key difference between clustering algorithms based on a "uncorrected" distance matrix and phylogeny is the latter is based on an explicit model to accommodate reversion mutations. The real problem is that there are only 4 bases, so randomly there is 1/4 chance that a mutation at a given position will mutate back to the original, e.g. A->C->A. Observed mutational differences = 0, expected (real) mutational differences = 2. What phylogeny is concerned about is reconstructing that "2". The problem is significant because almost every gene has regions of rapid mutations and regions of low mutations.
The primary way of doing this is via the Jukes-Cantor correction and is foundational in all nucleotide trees expect "p-distances". If nucleotide divergence is lower than 75% then it is possible to estimate the expected number of mutations using the observed via the JC correction. Moreover, when combined with a method of estimating rate variation within a gene, (usually the discrete gamma distribution), JC correction is very effective at recovering the "true tree". This is because fast evolving homologous sites are grouped together - big correction via JC, slow evolving sites are grouped together - small correction via JC. Other approaches to improving the JC correction are via identifying the bias between purine to purine mutations and pyrimidine to pyrimidine and purine to pyrimidine mutations.
The importance of JC was demonstrated by simulation studies of 4 taxa ((a,b), (c,d)), if two strains evolved very quickly when there sister lineages evolve slowly a clustering algorithm will report the fast strains are a sister group, i.e. ((b, c),(a, d)). If the JC method is implemented via maximum likelihood (or Bayesian) then it correctly recovers the true tree ((a,b), (c,d)). The artefact is known as long branch attraction.
Clustering algorithms based on a distance matrix that implements the JC correction, tend to perform poorly for long branch-attraction artefacts. This doesn't mean the correction is useless, just not particularly powerful. The issue is that distance matrix methods don't "stick" to the model and the clustering will introduce a layer of inaccuracy. Normally presentation of a "corrected" distance matrix combined with a clustering algorithm will require bootstrapping (resampling with replacement) to assess whether a given cluster is supported. Parameterised distance matrices, using neighbour-joining clustering in in-conjunction with a bootstrap are considered okay.
The @user172818 mentioned parsimony and this method is considered less reliable because it can't implement a JC correction. IMO it is possible weighted parsimony could make a "come-back", but it would be really complicated to implement a biological weighting method and would require extensive, independent calculations.

A great question, though a little ambiguous. I don't know what "general clustering algorithms" refer to. For biological sequences, building a tree can be thought as a way of clustering. Anyway...
There are different tree building algorithms. Max parsimony (MP), max likelihood (ML) and bayesian algorithms directly take sequences as input. They are distinct from distance based clustering.
Then there is a class of distance based algorithms in phylogenetics. They start from an all-pair distance matrix and aim to find a tree that is most compatible with the matrix. Among them, UPGMA is largely hierarchical clustering. Neighbor-joining is somewhat similar to hierarchical clustering, but it builds unrooted trees. FastME uses a very different approach. Briefly, for each topology it tries to find the best branch lengths by weighted least square (for details, see its paper); the best topology minimizes total branch lengths. FastME finds the best topology via nearest neighbor interchange (NNI), closer to some ML algorithms (e.g. PhyML).

does it brings anything new at the level of clustering algorithm?

Personally, I think the distance based tree building methods are more sophisticated than naive hierarchical clustering, provided that the real relationship follows a tree topology.


text:
Is there a way to assemble contigs starting from a specific sequence?
My work involves searching for marker genes/fragments in metagenomic databases (like the Sequence Read Archive). Once I find these sequences, I would like to know more about the neighboring genomic region.
Is there a way I could assemble only sequences that create a contig which contains my region of interest? Contigs which don't contain this region are not useful to me. My organism of interest might represent a minority of the metagenome, and assembling everything in the dataset would use a lot of computing power.

Good question! I know that some of the assemblers for mtDNA are working on this principle - searching for reads mapping to a conserved mt gene and then extending the sequence, maybe you could tweak one of those tools.

This is one of the primary use cases for which spacegraphcats (preprint and code) was designed. The "neighborhood queries" discussed in the paper sound particularly relevant. I don't have any personal experience with spacegraphcats, but the run guide provides some examples of how to index the complete data set and how to query sequences of interest.

Yes, we'd be happy to chat (I'm one of the spacegraphcats authors). Just post questions over on our github repo!

You might also be interested in citations 21 and 22 in the spacegraphcats paper - metacherchant and a paper that does something similar.


text:
Is it advisable to remove X and Y chromosome genes in a bulk RNA-seq dataset at the level of the count matrix?
From this link remove X and Y chromosome genes in RNA-seq data using DESeq2 pipeline, I have learned that depending on context, it is perfectly valid to remove X and Y chromosomal genes in an RNA-seq dataset before doing differential expression analysis. However I only have access to the count matrix and not the bam files of the RNA seq data I am analysing. Is it advisable to do the X and Y chromosome gene removal at the level of the count matrix? I am talking about removing all rows that contain X and Y chromosomal genes from the .csv count matrix file read in to RStudio before proceeding to create the DESeq2 object.
Thank you in advance for your kind response.

First off, it can be valid, but why do you want to do so? In my lab we use them as a validation (does sample X that is a male, has gene expression of Y?). I think that removing them leads to more work with unclear benefit: what if there is a relevant gene in the X chromosome?

Why do you want to do this? It is almost never a good idea and will most likely not have the effect that you think it will (e.g., it won't remove gender differences).

@Devon Ryan :  Thank you very much for your response. Do you please have a supporting publication?

No one would ever bother writing a paper about something that is (A) almost never done and (B) rarely useful.


text:
VEP output SIFT_score unclear
We have been experimenting with VEP (Variant Effect Predictor). One of the meta data attributes that we are interested in is the SIFT score, indeed when we apply the dbNSFP plug we get a column containing the scores (named SIFT_score). However, I don't understand why there are sometimes dots or multiple values in the fields. For example, the gene ENSG00000196924 below has 5 transcripts:

The SIFT_score column contains several values, not 1 per transcript/rs-number...
Here is another example that confuses me (I added the SIFT_pred column this time):

There are two mutations, the lower one can be expressed in 4 transcripts and I thus understand that there can be 4 SIFT scores, but why are all for in given in every row? Is the first one the SIFT_score for the first transcript?
One last example, again 4 transcripts, but now 2 of the scores are dots, what does that mean?

I have been looking for quite some time now how to interpret this data, any help is appreciated.

The dbNSFP plugin from VEP accesses tables of data for each variant from dbNSFP and pulls out the values. dbNSFP provide their SIFT scores in that format: a score for every transcript affected by the variant, all on one line. The lookup is just for the variant, not for the variant/transcript combo, so they provide scores for all variant/transcript combos. You can also get a column that gives you a list of the transcripts or proteins (Ensembl_transcriptid or Ensembl_proteinid) in order so you know which score goes with which transcript.
A better way to get SIFT scores with VEP is to get them directly from VEP, rather than using dbNSFP. This will get the SIFT score that goes with the transcript on the line with the relevant transcript.

I am guessing the dots are there for cases where dbNSFP doesn't have a value for the relevant transcript, right?

Yes, that's it. Could be that the variant isn't missense in that transcript.

The first gene you mention, ENSG00000196924, actually has 6 transcripts (link to the VarSome.com page of variant rs371839875), not 5. It's just that one of them is non-coding:

So the Sift scores you see are indeed one per transcript, it's just that there are 6 because dbNSFP also includes a score for the non-coding transcript of the gene.
The dots are just there as placeholders, they mean there was no value associated with that transcript. Many tools will show some sort of symbol instead of an empty field both for clarity and for practical technical reasons.
Visiting the variant's page on VarSome gives you a clearer picture since we collapse the identical scores and also include the converted rankscore provided by dbNSFP so you can have a single number for your variant:

Disclaimer: I work for the company behind VarSome, but it's a free tool. You need to pay to annotate VCF files (unlike the 100% free VEP), but it's free to use as a lookup tool for single variants.


text:
How to write the subclusters in file?
I have constructed a gene co-expression network from RNA-seq data. The network file is in edge list format of memory around 1gb which was created by calculating Pearson correlation of each gene pairs and the gene pairs which having correlation >95% were selected to create the edge list.
I have clustered this gene network (edge list) using igraph R package "cluster_louvian" community detection algorithm and obtained 534 subclusters.
c <- cluster_louvian(g)
c$membership

to get largest subcluster
sub1 <- induced.subgraph(g, membership(cl) ==order(-sizes(cl))[1])

I am finding difficulty to write the subcluster into a text file as an edge list
How can I save each subcluster in separate files (text files) as edge list?

Hi, could you explain better what software are you using? (Where is that cluster_louvian community detection?) Also in which format do you have these data? What have you tried to separate the subclusters? In what format do you want that edge list?

Now is more complete. However what have you tried? Could you give us some example data ?


text:
What does "fetching by region is not available for SAM files" mean?
I am used to gzip/biopython solutions when dealing with sequencing data, but now I wish to switch to more elegant pysam. So I looked at the manual, but ran into quite bizarre troubles with the first couple of lines using my bam file
import pysam
samfile = pysam.AlignmentFile("3_Tms_1_mapped.bam", "rb")
for read in samfile.fetch('3_Tms_b3v08_scaf000159'):
    print(read)
samfile.close()

returns ValueError: fetching by region is not available for SAM files. Well, the file is bam. I tried to google the error but the only hits I found are the lines in the source code of pysam that check if the file is bam/cram or sam, so somehow pysam thinks that my bam is a sam. How can I convince it otherwise? I have also noticed that the manual is for python 2.7, that's maybe where the problem comes from...

Your 3_Tms_1_mapped.bam file, despite its filename extension, is in fact a bgzipped SAM file. You can verify this using htsfile, which is a small utility packaged with HTSlib:
$ htsfile 3_Tms_1_mapped.bam
3_Tms_1_mapped.bam: SAM version 1.3 BGZF-compressed sequence data

(For files that really are in BAM format, it reports BAM version 1 compressed sequence data.)
So the error message is accurate in this case.

Hmm, does it mean that this link is wrong? http://biopython.org/DIST/docs/api/Bio.bgzf-module.html I sort of thought that it's the same thing...

No. (Was there a particular part of that you think might be wrong?) But there is a difference between BGZF-compressing the plain text format, and a BAM file (whose decompressed underlying stream is a bespoke binary format that is different from the plain SAM text).

I think I get it now. I got confused by the aims to this module, they write "In addition to being required for random access to and writing of BAM files", so I got an impression that I can simply use this module to write BAM files. This was actually really helpful answer, thanks!

@KamilSJaron you might want to consider accepting this answer instead of mine since this was the real issue.

That isn't actually a bam file as John Marshall figured out. I am keeping the rest of my answer since it could be useful to someone else, but the issue here was that you had a compressed (bgzipped) sam file and not an actual bam file and that's why you were getting that error. When I sorted your file in preparation for indexing it, I converted to a bam which is why the rest of my answer worked.

You don't have the index file for the bam file you're using. I used this script on the file you linked to (changing the name so that it corresponds to the right file and a contig in that file):
#!/usr/bin/env python3

import pysam
samfile = pysam.AlignmentFile("3_Tce_1_mapped.bam", "rb")
for read in samfile.fetch('3_Tce_b3v08_scaf005149'):
    print(read)
samfile.close()

The directory I ran it in had:
$ ls 3*
3_Tce_1_mapped.bam

And I got the error you described:
$ foo.py
Traceback (most recent call last):
  File "/home/terdon/scripts/foo.py", line 5, in <module>
    for read in samfile.fetch('3_Tce_b3v08_scaf005149'):
  File "pysam/libcalignmentfile.pyx", line 1107, in pysam.libcalignmentfile.AlignmentFile.fetch
ValueError: fetching by region is not available for SAM files

However, indexing the bam file fixed it:
$ samtools sort 3_Tce_1_mapped.bam > 3_Tce_1_mapped.sorted.bam
$ mv 3_Tce_1_mapped.sorted.bam 3_Tce_1_mapped.bam
$ samtools index 3_Tce_1_mapped.bam
$ ls 3*
3_Tce_1_mapped.bam  3_Tce_1_mapped.bam.bai

$ foo.py | wc
227    2724   16725

So just sort and index your files before attempting to seek in them. Which makes sense since the index's job is primarily to allow seeking.


text:
A GTF annotation for GRCH375d version
I have been given a bunch of .bam files, that come from the alignment on GRCh37d5 of human genome by Star tool. I want to extract raw read counts from these bam files by FEATURECOUNTS tool so I need a GTF. I googled a lot but I did not find GTF for GRCh37d5. I used another GTFs and I obtained weird features. Please one of you help me in finding this specific GTF. I converted a bam to fastq and I used Salmon to quantify gene expression by a transcripome.fasta but rate of mapping was 0.06

It's helpful to know that GRCh37 refers to the Genome Reference Consortium assembly of the human genome, version 37 (should be equivalent to hg19 from UCSC). I've never seen "5d" though. Often times GRC will introduce "patch" updates with minor changes that don't change the primary coordinate system—that's what p13 refers to in kayvee's answer below.

@DanielStandage Might be of interest to you: https://www.biostars.org/p/328824/#330791

The most recent annotation of GRCh37 assembly by NCBI RefSeq is here: https://ftp.ncbi.nlm.nih.gov/genomes/Homo_sapiens/GRCh37.p13_interim_annotation/
You can download the annotation in GFF3 format from that path and use it with featurecounts tool by first converting it to the SAF format using the commands shown below:
$ zgrep -v '^#' annotation.gff.gz \
    | awk 'BEGIN{FS="\t";OFS="\t"}($3=="exon"){print $1,$4,$5,$7,$9}' \
    | sed -r 's/ID.*GeneID:([0-9]*).*/\1/g' \
    | awk 'BEGIN{FS="\t";OFS="\t"}{print $5,$1,$2,$3,$4}'
    > annotation.saf

You should confirm that the chromosome sequence identifiers in the BAM file and the SAF file match. NCBI annotation has identifiers in the accession.version format whereas UCSC uses chr# format and Ensembl uses #. For this, you can use the -A parameter from featureCounts. Assuming that your BAM file uses UCSC-style identifiers, you can generate a 'chromosome name alias file' that can be used with featureCounts as follows:
## get NCBI assembly report for GRCh37
wget -nv ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.24_GRCh37.p12/GCF_000001405.24_GRCh37.p12_assembly_report.txt

## generate chr_alias.csv file
grep -v '^#' GCF_000001405.24_GRCh37.p12_assembly_report.txt \
    | dos2unix \
    | awk 'BEGIN{FS="\t";OFS=","}($10!="na"){print $10,$7}' \
    > chr_alias.csv

I didn't know about "SAF" format. I use gff with FeatureCounts. For the record, according to https://www.rdocumentation.org/packages/Rsubread/versions/1.22.2/topics/featureCounts: "The SAF annotation format has five required columns, including GeneID, Chr, Start, End and Strand. These columns can be in any order".


text:
Consensus sequence from SAM or BAM file?
I am trying to perform reference-based assembly. Most of the tutorials teach how to create a bam file and view alignemnts in IGV or Tablet. But, I want a assembled genome sequence in fasta format. How can I do that ? Is there any tool or linux syntax available for this purpose ?

Cross-posted on biostars: https://www.biostars.org/p/374352/

A few years ago I used AlignGraph (paper and code) for a task like this. I gave AlignGraph my raw reads, so I'm not sure if it accepts SAM/BAM input or if you have to go back to the original Fastq reads.
I'm not sure if any related/improved tools have been published since AlignGraph. Internet searches for "reference-guided assembly" return mostly transcriptome assembly tools. It might be worth looking for articles that cite the AlignGraph paper.

It seems like AlignGraph takes reads in fasta format. Does it accept fastq format ? What tool do you use to convert fastq to fasta format ?

cat reads.fastq | paste - - - - | cut -f 1-2 | sed 's/^@/>/' | tr '\t' '\n' > reads.fasta. You can replace cat with gunzip -c if your Fastqs are compressed.


text:
Post docking evaluation
What to do after protein-ligand docking? Is there any necessary or enhancement step recommended after docking to validate the confirmation of highest score is the best? For your information, I am using CCDC GOLD for docking and GOLD score as my scoring function.

Why do you need to to a docking? Do you have a biological question? For example, I want to know if this drug docks to this receptor. To validate the docking you usually test it in real life...

If you really have to, you can try different docking algorithms and scoring FF. However, you might get subtle differences due to how your ligand is parameterised.

just a follow up to the answer, what if we don't want to proceed with the assays and just want to stick with the bioinformatics. What else can I do to make my filtered drug molecule a stronger candidate?

You can run a Molecular Dynamics simulation to see if the pose is stable and whether the protein changes/adapts it’s conformation.
You can calculate the ligand efficiency, the docking score divided by the number if heavy atoms, and rank the compounds accordingly.
You can also compare the scores of known ligands to the experimental values and see if the docking the scores are predicting known ligands or correlate with the experimental ones.
At some point you have to make a selection and purchase or synthesize the compounds. How many compounds can you test experimentally? How many can you purchase? What is the purpose of these ligands going to be? Drugs? Did you screen the correct set of ligands e.g. a library of lead compounds that are ready to purchase for cheap, that would be a good start.

To further this answer, a form of steered MD that I have found promising is (dynamic undocking)[https://www.nature.com/articles/nchem.2660]: your ligand gets yanked out. The output is not super clear though and requires some manual checking.


text:
Parsing .vcf file for this information
I have a .vcf file
https://www.dropbox.com/s/8v73nppwg3a1tnd/LP2000109-DNA_A01_vs_LP2000103-DNA_A01.SVannotated.txt?dl=0
with this header
##startTime=Fri Mar 29 16:46:32 2019
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  NORMAL  TUMOR
1   54586   .   T   C   .   PASS    DP=39;MQ=50.55;MQ0=0;NT=ref;QSS=48;QSS_NT=48;ReadPosRankSum=1.92;SGT=TT->CT;SNVSB=0.00;SOMATIC;SomaticEVS=10.83;TQSS=1;TQSS_NT=1    AU:CU:DP:FDP:GU:SDP:SUBDP:TU    0,0:0,0:20:0:0,0:0:0:20,20  0,0:6,6:18:0:0,0:0:0:12,13
1   103241  .   C   T   .   PASS    DP=120;MQ=24.94;MQ0=35;NT=ref;QSS=47;QSS_NT=47;ReadPosRankSum=2.09;SGT=CC->CT;SNVSB=0.00;SOMATIC;SomaticEVS=9.44;TQSS=2;TQSS_NT=2   AU:CU:DP:FDP:GU:SDP:SUBDP:TU    0,1:32,47:33:1:0,0:0:0:0,5  0,

The "DP" field in the vcf shows the depth of the individual samples; So in this file, the first locus has the following format fields:
AU:CU:DP:FDP:GU:SDP:SUBDP:TU    0,0:0,0:20:0:0,0:0:0:20,20  0,0:6,6:18:0:0,0:0:0:12,13

So according to this (DP field of normal and tumor samples), normal sample has a depth of 20 and tumor sample has a depth of 18.
So how I could extract the read depth for all loci as described for the first position? The desired output would be like this [note that the VCF is taken from my own data but the table is my desired format that I don't know how to get that from my own data. chr have been added manually because my reference genome is hg19]:
  Sample   Type   CHROM       POS    REF  ALT       Tumor_Depth Normal_Depth
  CHC2432T SNV    chr1  102961055    G    A                  64           62
  CHC2432T SNV    chr1  105492588    A    T                  66           73
  CHC2432T SNV    chr1  108628724    C    T                  45           54
  CHC2432T SNV    chr1  109692113    G    T                  53           29
  CHC2432T SNV    chr1  109692114    G    T                  53           31
  CHC2432T SNV    chr1  120676701    T    C                  48           87

Your vcf file is truncated (see tail -n1 _double_filtered_013_pre_indels.vcf). I assume that's just the copy you uploaded to dropbox, but check your actual file to be sure.

To extract the DP fields from a VCF file, you could use a tool like bcftools query:

Extracts fields from VCF or BCF files and outputs them in user-defined format.

You could start from something like this:
bcftools query -Hf 'CHC2432T\t%TYPE\t%CHROM\t%POS\t%REF\t%ALT[\t%DP]\n' file.vcf

Sorry, how I can add RU=TAC part to tumour sample?

You can do the extraction part with the GATK tool VariantsToTable, as described here:
https://gatk.broadinstitute.org/hc/en-us/articles/360041414592-VariantsToTable
The usage example from that doc:
gatk VariantsToTable \
    -V input.vcf \
    -F CHROM -F POS -F TYPE -GF AD \
    -O output.table

would produce a file that looks like:
 CHROM  POS        TYPE   HSCX1010N.AD  HSCX1010T.AD
 1      31782997   SNP    77,0          53,4
 1      40125052   SNP    97,0          92,7
 1      65068538   SNP    49,0          35,4
 1      111146235  SNP    69,1          77,4

So you might still need to reorder columns etc but that should allow you to at least get the values out in a tabular format that will be easier to work with.