text:
How can the cell line contribution be estimated from RNASeq data?
Using a laser-capture microdissection of cells a group of cells stained with the marker of interest was sequenced. In another cohort of patients (this is all  human liver tissue) the whole tissue was sequenced (RNA-seq in both cases)
Can I estimate the contribution of the cells marked in the whole liver ("weight of these cells" in the liver in words of my PI)?
My gut feeling is that it can't be done this way, it would require both single cell sequencing and whole tissue sequencing to estimate the contribution of each cell line. But perhaps there is a tool that given the cell lines or the main expression of other cells lines it can be compared to using GSVA or some similar tool.

Have you looked at tools for admixture estimation (that's what you're doing)? How precise do you need the estimate to be?

I haven't heard of admixture estimation, so I haven't look for tools with that keyword. I don't have a requirement of precision, the more, the better :D. But I suspect my data is not too much good (I have just 6 technical replicates of the laser microdissection) so I can't expect much.

There are couple of computational methods which try to do this (I never used them, so no experience):

CellMix, based on sets of marker gene lists
Subset Prediction from Enrichment Correlation, which is based on correlations with subset-specific genes across a set of samples.
Cell type enrichment, which uses our highly expressed, cell specific gene database
Cell type-specific significance analysis using differential gene expression for each cell type

You might have to get some reference expression levels from public databases or papers for some of the methods.
One thing to keep in mind: you cannot really compute the cells proportion, only RNA proportion. If you have a good reason to assume that RNA quantity
per cell is very similar, this is a good proxy for the cells proportion in a tissue.

Nice references. I'll have that limitation in mind. Thanks

Cell deconvolution is mentioned in this Biostars post, which mentions CIBERSORT for immune cell mixes, and the Bioconductor package DeconRNASeq.
As far as I'm aware, it is only possible at best to get proportional representation for transcript expression from standard high-throughput sequencing results, because the sequencers and sample preparation workflow are designed in such a way that the same number of reads are output regardless of the input amount.

CIBERSORT sounds like a nice tool, worth a try.

You can download public single-cell RNA-seq data or bulk RNA-seq for purified primary cells for reference and then use CIBERSORT or MuSiC to deconvolve.
One thing to notice: to make bulk RNA-seq and single-cell reference compatible, TPM from RNA-seq can be compared with CPM in 3' (10X,Drop-seq etc.) single-cell data and TPM in full-length single-cell data (SMARTseq etc.)

Actually a month ago, I was searching for the same question and I found a few packages that can be used. First of all the functionality of them can be cell proportion estimation, cell specific gene expression estimation, or both of these functions. Secondly, their primary input, as expected, is bulk RNA-seq data and some also want a single cell RNA-seq data as secondary input.
The most popular one is the child of famous cibersort, CIBERSORTX. I went through their website and read the tutorials. They have several example datasets that can be used for learning their platform. However, the two biggest flaws were that they haven't shared their source code which means we can't reproduce the result on our own systems and we have to upload all of our data to Stanford servers (security concerns!). The second problem is that there is no API access to their server and everything should be done online which means very low compatibility with other tools.
The idea, using deconvolution and deep learning to infer subpopulation of cells from bulk RNA-seq, is great. So, I searched to see if others have also done something similar and maybe even better. I found three other tools:

MuSiC:   https://www.nature.com/articles/s41467-018-08023-x
Published in nature communication and a few months older than cibersortx. I haven't studied their algorithm yet and so I can't say if their method is better than cibersortx or not. However, the part that makes the difference is that all of their works, including their source code, is open-source. So, I can access their code and we don't need to upload our unpublished data into some not safe online server (unlike cibersortx). The second good thing is that it is implemented in R language which I use for most of my RNA-seq data analysis, so it is easily compatible with other tools.

Scaden : https://advances.sciencemag.org/content/6/30/eaba2619.full
This one is also open-source and accessible easily. It is written in Python and runs in the shell. So, it is somewhat compatible with other tools (not as compatible as MuSic). We again can process our data locally and don't need to upload our data to any servers (which is good). They claim that their method has a better correlation with the real scRNA-seq data than both cibersortx and MuSic (not always!!). Running their tools seems easy based on their tutorials (I haven't run it on my system yet!). However, their tutorials seem to be a little confusing and it seems like it is still in pre-release version and they have bugs!!

CDSeq:
https://github.com/kkang7/CDSeq_R_Package
Finally, I found the CDSeq package which doesn't need any scRNA-seq input to find single cell proportions and cell-specific gene expression in bulk data. It just gets your buld counts as input and info about very few parameters and does the whole thing automatically, and suprisingly good, I checked! If you could find gene expression pattern of your specific cells and also use it as an input, it is optional, it will give you even a better estimation.


text:
Quantifying reads mapping to multiple loci
I have been using STAR for our RNA-Seq samples. The final.out log file reports percentage of uniquely mapped reads along with percentage of reads that map to multiple loci (less than or equal to 10) and percentage of reads mapping to too many loci (greater than 10). However, I want to break down the multiple loci part to individual counts: Reads mapping to 2 locations, 3 locations, 4 locations .. 10 locations.
The NH tag seems to be used by STAR. However a naive read counting approach results in it reporting more number of reads than total reads.
For example, my final.out looks like this:
   Mapping speed, Million of reads per hour |       1403.36

                      Number of input reads |       53015978
                  Average input read length |       26
                                UNIQUE READS:
               Uniquely mapped reads number |       368916
                    Uniquely mapped reads % |       0.70%
                      Average mapped length |       26.45
                   Number of splices: Total |       1057
        Number of splices: Annotated (sjdb) |       0
                   Number of splices: GT/AG |       802
                   Number of splices: GC/AG |       1
                   Number of splices: AT/AC |       0
           Number of splices: Non-canonical |       254
                  Mismatch rate per base, % |       0.31%
                     Deletion rate per base |       0.00%
                    Deletion average length |       1.45
                    Insertion rate per base |       0.00%
                   Insertion average length |       1.00
                         MULTI-MAPPING READS:
    Number of reads mapped to multiple loci |       45766732
         % of reads mapped to multiple loci |       86.33%
    Number of reads mapped to too many loci |       3757890
         % of reads mapped to too many loci |       7.09%
                              UNMAPPED READS:
   % of reads unmapped: too many mismatches |       0.00%
             % of reads unmapped: too short |       5.89%
                 % of reads unmapped: other |       0.00%

Counting histogram of number of positions a read maps to using pysam:
def get_reads_hist(bam):
    bam = pysam.AlignmentFile(bam, 'rb')
    counts = Counter()
    for query in bam.fetch():
        nh_count = Counter(dict(query.get_tags())['NH'])
        counts += nh_count
    return counts

results in
Counter({1: 330606,
     2: 86772164,
     3: 329,
     4: 38083,
     5: 31,
     6: 1094,
     7: 129,
     8: 50,
     10: 50})

The count 1 reads are fine even though they do not match the counts in final.out file since I am counting a certain category of reads (say those mapping to tRNA only), but the reads mapping to 2 locations are highly overestimated. Why is that?

You almost had the correct python code already, you just need to filter out secondary alignments:
def get_reads_hist(bam):
    bam = pysam.AlignmentFile(bam, 'rb')
    counts = Counter()
    for query in bam.fetch():
        if query.is_secondary:
            continue
        nh_count = Counter(dict(query.get_tags())['NH'])
        counts += nh_count
    return counts

source: Stack Exchange

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

sort: votes


text:
How to manipulate protein interaction network from String database in R?
How can I manipulate protein-interaction network graph from the String database using STRINGdb bioconductor package and R?
I have downloaded the whole graph for Homo sapiens from STRING, which has about 20.000 proteins.

How do I read the file using that package?
How do I filter things I don't need? Supposing that I want to keep tumor data, as an example.

Cross-posted on biostars...again

I don't think it is a problem to cross-post. It would be good to post a link to the best answer on the other site when such an answer is defined.

@bli that what I'm going to do once an answer is found

That's the question how do I do to filter things I don't need I already downloaded the whole graph for homosapiens from the download section in STRING .

I think the easiest way is to download the graph using STRINGdb.
library(STRINGdb)
string_db <- STRINGdb$new(version="10", species=9606,
                          score_threshold=400, input_directory="" )
full.graph <- string_db$get_graph()

Now you can use igraph, to manipulate the graph. Let's assume you want to take 200 proteins with the highest degree, i.e. number of edges they have.
library(igraph)

# see how many proteins do you have
vcount(full.graph)

# find top 200 proteins with the highest degree
top.degree.verticies <- names(tail(sort(degree(full.graph)), 200))

# extract the relevant subgraph
top.subgraph <- induced_subgraph(full.graph, top.degree.verticies)

# count the number of proteins in it
vcount(top.subgraph)

How to get disease specific genes?
There's no GO annotation for cancer or Alzheimer's disease. It is out of scope of the GO consortium.
What you can do, you can either take KEGG Pathways annotation, or manually select list of relevant GO-terms. Or acquire the list from one of the papers. For example annotation term 05200 corresponds to the cancer KEGG pathway. You can easily retrieve proteins associated with the annotation:
cancer.pathway.proteins <-
    string_db$get_term_proteins('05200')$STRING_id

And then perform subgraphing as described above.
Alternatively you can try to get an enrichment score for an every gene given it's neighbors (the way enrichment is shown on the string-db website). Then you can keep only those having top enrichment scores. Probably get_ppi_enrichment_full or get_ppi_enrichment functions will help you to do that.

If we want to keep tumor , cancer , alzheimer data , I mean filtering by disease. can we do that ? knowing that we want to keep the maximum of nodes

@I'm getting errors Error in induced_subgraph(graph, top.degree.verticies) : Not a graph object

That's probably because graph is a function in igraph package. Make sure you're using a correct variable name.

is there a way to export the graph is csv ?

Please, check the very good documentation of the igraph package. You can export graph using write_graph. Also, as you can notice STRINGdb downloads the full graph to your computer, therefore you probably can download it directly from the database website.

Let us continue this discussion in chat.

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


text:
Is there any publicly available multi-omics dataset?
I am searching for a multi-omics dataset which may include genomics, transcriptomics and proteomics (eg. snyderome).
I need such a dataset for an introductory data-exploration purpose. So it can be from any host organism and any experimental methods/conditions.
Is there any database/repository where I can find it?

Do you really need all 3 or will any 2 of those work? For "genomics", do you want histone marks (or methylation) or just DNAseq (e.g., for variant calling)?

genomics (DNAseq) would be essential. rest any layers (at least two) would work.

Please [edit] your question and clarify what you need. Does everything need to be from the same cell or individual? Are you looking for the set of expressed proteins, transcribed RNA and the associated genome sequence? Or maybe those of a specific condition? What's parallel genomics?

I wanted to make an suggestion for rephrasing the question, but apparently I made the edit. If it does not say now what you wanted, feel free to change it back...

There is a database called OmicsDI, where one can search for multi-omics datasets.
Here's a link of the associated publication (Perez-Riverol, Yasset, et al. "Omics Discovery Index-Discovering and Linking Public Omics Datasets." bioRxiv (2016): 049205.) for more details.

There are several datasets available on GEO, though you do have to search for them. For example, here are three data sets that have both Illumina methylation and gene expression microarray profiling:

GSE87650
GSE72874
GSE56047

Next to OmicsDI the EBI has a special repository for multi-omics datasets: https://www.ebi.ac.uk/biosamples
It links the different datasets between repositories, ie. PRIDE for MS/MS based data and ArrayExpress for RNASeq data.


text:
Remapping genomic coordinates to account for indels
I'm interested in obtaining coding sequences of my favourite gene in all individuals from the 1000Genomes (and similar projects). I use GATK to get the right subset of variants, vcf-consensus to map these variants onto the reference genome and finally samtools to extract the individual exons. This works fine if the variants are SNPs but if there are any indels, this changes the coordinates of exons and I end up getting the wrong region. Is there any generic way of remapping genomic coordinates to account for the changes created by indels?

I think that you need a LiftOver Chain file to transform your coordinates. You can obtain such a file using bcftools consensus with the -c parameter:
-c, --chain <file>         write a chain file for liftover

Then you can use it to transform coordinates in various genomic formats using CrossMap.

Nice, +1. It would be great if you could expand this to include an example command the OP could run, the steps needed etc. The Stack Exchange sites have very strict guidelines about what is an answer (that's why I've been bugging you so) and ideally, we want that an answer provide all the needed information to solve the question asked. In other words, please don't post pointers (not even very helpful ones like this) as answers. Either just post a comment, or flesh the pointer into a full answer.

Maybe is just me, but his answer seems fine, even to SE standard, but it's just me.

@nuin agreed, which is why I upvoted. I just think it would be even better with more detail, that's all.

No hard feelings.

If the goal is to extract the consensus sequence for given regions, this is how it works in 2019.
bgzip and index your vcf file.
$ bgzip -c input.vcf > input.vcf.gz
$ tabix input.vcf.gz

Create a regions.txt which contains one region per line in the format chr:from-to. If you have already a bed file, you can use this little awkscript to create it:
$ awk '{print $1":"$2+1"-"$3}' input.bed > regions.txt

Now run this combination of samtools and bcftools:
$ samtools faidx -r regions.txt genome.fa|bcftools consensus input.vcf.gz -o consensus.fa

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


text:
Have DNA motifs 6-12bp long, trying to get conservation scores
I have about 200 short nucleotide motifs (6-12 bp in length) from the human genome, and I'm trying to see how conserved they are across vertebrates.
I was thinking that I'd need to make a bed file for each motif that lists all of its occurrences in the human genome. From there, I could map the beds to a bigwig files of PhastCons scores (essentially doing the reverse of what the PhastCons software was designed to do). Does that sound like the best approach?
I'm getting stuck at the step of going from motifs to bed files. I've tried using BLAST to find all occurrences of motifs, but their short length is causing issues.
I've tried messing with the e-value threshold, word size, and filter parameters, but I still don't get any hits.
Is there a work-around for this issue, or should I just rethink my entire approach?

I did something similar, but with a list of well defined n-mers instead of "motifs", all with the same value of n. In case this can be of any help, the code I used to make the bed file is the following: https://bitbucket.org/blaiseli/conservation_mirna_targets/src/20b7a0894b957e69912f637ef2d0493779036029/trie.py?at=master&fileviewer=file-view-default
The code might not be super memory efficient.

In case you have only ACGT in your motifs
The short motifs make it sound as if you are in the business of looking for a kmer counter. You can choose to use existing software or build your own.

Using existing software might be your easiest path. An older post from 2014 will probably give you a first idea what's out there: http://homolog.us/blogs/blog/2014/04/07/kmer-counting-a-2014-recap/ . Note that a couple of algorithms mentioned there got successors, so it is worthwhile digging a bit around. The small kmer size will make most of them usable for your needs.
As the maximum size of your kmers is comparatively small (12 nt need 24 bits, i.e., max 16.7 million entries in your kmer table), you should be able to easily roll your own kmer counting in about any language you like and on about any of nowadays computer. The pseudocode section on the Wikipedia entry for kmers will give you first pointers for that. Might be a bit more work, but maybe more flexible depending on your needs.

In case you have IUPAC bases (N, W, etc.) in your motifs
I don't know any any pre-existing software doing what you need. I could imagine that the short motifs make using regular expressions doable for this kind of search, but I may be wrong. Testing this should be easy in a simple script as all major programming languages have modules or libraries for REs. Even if it should take a couple of hours to run on your data set, that would be good enough for a one-off calculation.

To scan motifs in a genome (or database) I would use FIMO which will give you the exact locations of these motifs in your genome.
Once you have the locations, you can use a phastCons bigiwig from UCSC to calculate the basewise conservation scores. However, please remember that phastCons scores are smoothed across windows and it might not be the best metric if you are trying to compare the conservation levels at your motif matching sites as compared to the sequences flanking them.
I wrote a package a while back to do this, including doing de-novo motif discovery. However, it might be an overkill for your use case.

I've used FIMO before actually, but for comparing motifs against PWMs (from JASPAR). Might be a dumb question, but is it pretty easy to compare against a genome rather than a database of PWM's?

@EricBrenner You can simply upload a sequence file in the web version.

Where transcription factor binding sites are concerned, in our papers we have used per-base phyloP conservation data instead of smoothed phastCons scores.
We use BEDOPS bedmap to map scores over multiple, padded binding sites (BED-formatted) for a given motif model. The phyloP scores are WIG files obtained from UCSC goldenpath and converted to BED via wig2bed.
The resulting matrix of binding site intervals and their per-base scores can be turned into a ranked heatmap or aggregated to determine average per-base conservation for a motif model.
While TF binding sites innately show high information content and therefore high conservation, matrices can be further sorted by score-maps of ChIP-seq or DNaseI-seq tag density.
This can help filter for interesting low-noise, high-signal patterns and is useful when the matrix is windowed and the window includes low-information regions with a few high-information residues outside of the main motif — as is the case with CTCF, for instance.


text:
How does the BWA-MEM algorithm assign its mapping qualities?
Is there any resource (paper, blogpost, Github gist, etc.) describing the BWA-MEM algorithm for assigning mapping qualities? I vaguely remember that I have somewhere seen a formula for SE reads, which looked like
$C * (s_1 - s_2) / s_1,$
where $s_1$ and $s_2$ denoted the alignment scores of two best alignments and C was some constant.
I believe that a reimplementation of this algorithm in some scripting language could be very useful for the bioinfo community. For instance, I sometimes test various mapping methods and some of them tend to find good alignments, but fail in assigning appropriate qualities. Therefore, I would like to re-assign all the mapping qualities in a SAM file with the BWA-MEM algorithm.
Btw. This algorithm must already have been implemented outside BWA, see the BWA-MEM paper:

GEM does not compute mapping quality. Its
  mapping quality is estimated with a BWA-like algorithm with suboptimal
  alignments available.

Unfortunately, the BWA-MEM paper repo contains only the resulting .eval files.
Update: The question is not about the algorithm for computing alignment scores. Mapping qualities and alignment scores are two different things:

Alignment score quantifies the similarity between two sequences (e.g., a read and a reference sequence)
Mapping quality (MAQ) quantifies the probability that a read is aligned to a wrong position.

Even alignments with high scores can have a very low mapping quality.

Unfortunately I don’t know the answer for BWA-MEM (since it differs from BWA!) but pretty much all other tools are described here: https://sequencing.qcfail.com/articles/mapq-values-are-really-useful-but-their-implementation-is-a-mess/

Maybe this page can help http://genome.sph.umich.edu/wiki/Mapping_Quality_Scores . Once you know the best and alternative positions where a read can align (or even the best and second best only?) it's not too difficult to implement I guess.

Yes, there bwa-mem was published as a preprint

BWA-MEM’s seed extension differs from the standard seed extension in two aspects. Firstly, suppose at a certain extension step we come to reference
  position x with the best extension score achieved at query position y.
...
Secondly, while extending a seed, BWA-MEM tries to keep track of the
  best extension score reaching the end of the query sequence

And there is a description of the scoring algorithm directly in the source code of bwa-mem (lines 22 - 44), but maybe the only solution is really to go though the source code.

Thank you for your answer. However, the question is more about assigning mapping qualities. Even reads with a very high alignment score can have a mapping quality equal to zero.

Have you checked the source code? line 22 - 44.

@KamilSJaron Wow, that’s terribly hard to understand. The actual code unfortunately isn’t better. :-( At any rate, could you update your answer to include this more prominently?

@ KamilSJaron I did and it is still not completely clear for me even in the easier case of single-end reads.

Well I had no intention to actually explain the score (since I really do not know and the question explicitly asked about resources). I just knew about the preprint and I also got the idea to look at the source code, where I found those 22 lines of maths that seemed to explain the scoring.

The MAPQ score appears to be calculated in the mem_approx_mapq_se function. I tried porting it to Python and inserting comments using ChatGPT 4, though I haven't tested it to ensure the output is identical to the original.
import math

# Constants
MEM_MAPQ_COEFFICIENT = 30.0
MEM_MAPQ_MAXIMUM = 60

def calculate_mapping_quality(alignment_options, alignment_region):
    # Determine the suboptimal alignment score, which is the second-best SW alignment score if available,
    # otherwise calculated from minimum seed length and match score.
    suboptimal_alignment_score = (alignment_region.suboptimal_score if alignment_region.suboptimal_score
                                  else alignment_options.minimum_seed_length * alignment_options.match_score)

    # Update the suboptimal alignment score to the higher of the two: current suboptimal score or tandem hit score.
    suboptimal_alignment_score = (alignment_region.tandem_hit_score if alignment_region.tandem_hit_score > suboptimal_alignment_score
                                  else suboptimal_alignment_score)

    # Early exit with zero mapping quality if the suboptimal score is greater than or equal to the alignment score.
    if suboptimal_alignment_score >= alignment_region.best_alignment_score:
        return 0

    # Calculate the maximum alignment length between the query and reference sequences.
    alignment_length = max(alignment_region.query_end - alignment_region.query_start,
                           alignment_region.reference_end - alignment_region.reference_start)

    # Calculate the identity of the alignment.
    identity = (1.0 - (alignment_length * alignment_options.match_score - alignment_region.best_alignment_score)
                / (alignment_options.match_score + alignment_options.mismatch_penalty) / alignment_length)

    # Compute the mapping quality based on the alignment score, identity, and other factors.
    if alignment_region.best_alignment_score == 0:
        mapping_quality = 0
    elif alignment_options.mapQ_coef_length > 0:
        # Adjust mapping quality based on alignment length and identity.
        length_factor = (1.0 if alignment_length < alignment_options.mapQ_coef_length
                         else alignment_options.mapQ_coef_fac / math.log(alignment_length))
        length_factor *= identity * identity
        mapping_quality = int(6.02 * (alignment_region.best_alignment_score - suboptimal_alignment_score)
                              / alignment_options.match_score * length_factor * length_factor + 0.499)
    else:
        # Use the constant MEM_MAPQ_COEFFICIENT for mapping quality calculation.
        mapping_quality = int(MEM_MAPQ_COEFFICIENT * (1.0 - suboptimal_alignment_score / alignment_region.best_alignment_score)
                              * math.log(alignment_region.seed_coverage) + 0.499)
        mapping_quality = int(mapping_quality * identity * identity + 0.499) if identity < 0.95 else mapping_quality

    # Penalize the mapping quality for the presence of suboptimal hits.
    if alignment_region.suboptimal_hit_count > 0:
        mapping_quality -= int(4.343 * math.log(alignment_region.suboptimal_hit_count + 1) + 0.499)

    # Ensure the mapping quality does not exceed the maximum allowed value.
    mapping_quality = min(max(mapping_quality, 0), MEM_MAPQ_MAXIMUM)

    # Adjust the mapping quality based on the fraction of repetitive regions in the alignment.
    mapping_quality = int(mapping_quality * (1.0 - alignment_region.repetitive_region_fraction) + 0.499)

    # Return the final mapping quality value.
    return mapping_quality

```


text:
Is it possible to perform MinION sequencing offline?
I vaguely remember, that the original plan of Oxford Nanopore was to provide cheap sequencers (MinION), but charge for base-calling. For that reason the base-calling was performed in the cloud, and the plan was to make it commercial once the technology is established.
Of course, this limits the potential uses of MinION in the field, since huge areas do not have decent internet connection. Also, not all the data can be legally transferred to a third-party company in the clinical studies.
For example for the Ebola paper, they had to create a special version of their software:

An offline-capable version of MinKNOW, with internet ‘ping’ disabled
  and online updates disabled was made available to us by Oxford
  Nanopore Technologies specifically for the project

There are couple of third-party base-callers available today. I am aware of Nanocall and DeepNano, but since they are not official, it can be hard for them to keep-up with the latest versions of sequencers and cells.

Is it possible as of today to sequence offline without a special arrangement (like the Ebola one).
If not, what's the policy of Oxford Nanopore toward third-party base-callers? Are they going to help them, or to sue them eventually?

Short answer: yes, but you need to get permission (and modified software) from ONT before doing that.
... but that doesn't tell the whole story. This question has the potential to be very confusing, and that's through no fault of the questioner. The issue is that for the MinION, sequencing (or more specifically, generating the raw data in the form of an electrical signal trace) is distinct and separable from base calling. Many other sequencers also have distinct raw data and base-calling phases, but they're not democratised to the degree they are on the MinION.
The "sequencing" part of MinION sequencing is carried out by ONT software, namely MinKNOW. As explained to me during PoreCampAU 2017, when the MinION is initially plugged into a computer it is missing the firmware necessary to carry out the sequencing. The most recent version of this firmware is usually downloaded at the start of a sequencing run by sending a request to ONT servers. In the usual case, you can't do sequencing without being able to access those servers, and you can't do sequencing without ONT knowing about it. However, ONT acknowledge that there are people out there who won't have Internet access when sequencing (e.g. sequencing Ebola in Africa, or metagenomic sequencing in the middle of the ocean), and an email to <support@nanoporetech.com> with reasons is likely to result in a quick software fix to the local sequencing problem.
Once the raw signals are acquired, the "base-calling" part of MinION sequencing can be done anywhere. The ONT-maintained basecaller is Albacore, and this will get the first model updates whenever the sequencing technology is changed (which happens a lot). Albacore is a local basecaller which can be obtained from ONT by browsing through their community pages (available to anyone who has a MinION); ONT switched to only allowing people to do basecalling locally in about April 2017, after establishing that using AWS servers was just too expensive. Albacore is open source and free-as-in-beer, but has a restrictive licensing agreement which limits the distribution (and modification) of the program. However, Albacore is not the only available basecaller. ONT provide a FOSS basecaller called nanonet. It's a little bit behind Albacore on technology, but ONT have said that all useful Albacore changes will eventually propagate through to nanonet. There is another non-ONT basecaller that I'm aware of which uses a neural network for basecalling: deepnano. Other basecallers exist, each varying distances away technology-wise, and I expect that more will appear in the future as the technology stabilises and more change-resistant computer scientists get in on the act.
Edit: ONT has just pulled back the curtain on their basecalling software; all the repositories that I've looked at so far (except for the Cliveome) have been released under the Mozilla Public License (free and open source, with some conditions and limitations). Included in that software repository is Scrappie, which is their testing / bleeding-edge version of Albacore.

In short, yes offline basecalling is supported.
Oxford nanopore has made available an official offline basecaller, called albacore. If you have an account to the nanpore community website you can download it from here.
Concerning running it in the field, last time I checked live basecalling was still requiring a lot of computing power, so you might want to run the MinION in the field, but basecall once back a your lab / computing facility.
I'm not aware of any plans to sue the third-party basecallers, I think Oxford Nanopore is pretty positive about them :-)

Running the MinION in the field and base-calling back in the lab defeats the whole (well, part of) the purpose: namely, to get results instantaneously.

True in some cases, but sometimes you want to sequence samples in remote areas / places where it'd be impractical / impossible to go back to a lab with good quality DNA samples (i.e difficult access to dry ice for transport). The MinION makes it possible to sequence in the field, then analyse back in the lab.
EDIT: Also, for small sequencing run you can probably basecall with a good laptop (8+ cpus)


text:
How can I improve the yield of MinION sequencing runs?
This is a frequently-asked question within the nanopore community. Oxford Nanopore currently claims that they are able to generate run yields of 10-15 gigabases (e.g. see here and here), and yet it's more common to see users only managing in the 1-5 gigabase range.
So why the big difference in yield?

I attended a talk by Josh Quick at PoreCampAU 2017, in which he discussed some common barriers to getting both good sequencing yield and long read length. It mostly boils down to being more careful with the sample preparation. Bear in mind that the MinION will still sequence a dirty sample, it will just be at a reduced yield. Here are my notes from that talk:

The hardest thing about MinION sequencing is getting the sample in the first place
There are lots of different sample types and extraction methods
You can't get longer reads than what you put in; shit in leads to shit out
DNA is very stable when not moving, but very sensitive to lateral damage
The phenol chloroform method of DNA extraction is very good, and can be used with a phase-locked gel to make extraction easier
A simple salt + alcohol extraction might be the best method for extraction (because it involves the least amount of work on the DNA)
EDTA (e.g. as found in TE buffer, and many extraction kits) is not compatible with the rapid kit
The most consistently good Nanopore runs produced by Josh's lab were 1D ligations runs on R9.4; the best overall run was a phenol-chloroform extraction + rapid kit
John Tyson can tune himself out of a low-yield hole (via software)
Getting small numbers of short reads is very important
Suggested (and mostly untested) purification techniques: spin column (60-100kb); ethanol extraction (100-150kb), dialysis (150-250kb); low melting-point agarose plug (~1Mb, DNA extraction in situ)
The nanopore protocol input is in nanograms, but should really be stated as molarity; the kit expects about 0.2 pmol input
Picture molecules tethered to the surface of the membrane. You can then see that the flow cell density is independent of the sequence length
Tapestation, Qubit and Nanodrop are all a good idea; a DNA sample that can pass all three tests will work well: no short-length shoulders by Tapestation, sufficient DNA by Qubit, high purity by Nanodrop
RNA can interfere with sequencing; digesting RNA is recommended
Freezing DNA is a really bad idea. The ice crystals are very good at chopping DNA up into small pieces
DNA that is kept in the fridge is remarkably stable; can be kept for over two years (and probably indefinitely)

Update on yield:
We had a run in June 2018 which produced about 3 million reads from PCR-amplified mouse cDNA (read N50 of ~1kb), and that was in a situation where the it was the first sample prepared by someone, the input DNA amount was probably about 20% of what it should have been, and we suffered a dataloss incident. If that's what happened with a bad run, I have high hopes for our runs in the future.
Our next cDNA run in August 2018 produced about 7 million reads and a similar read N50 (i.e. total yield of about 7Gb). The run was interrupted by a Windows update during the first night (I had disabled automatic updates previously, but they had been re-enabled), and there have since been additional software and hardware updates since to improve sequencing yield. If we have a good flow cell and similar sample prep, I expect that the next cDNA run that we'll do should produce over 8 million reads, and possibly over 10 million.

My bet is that the nanopore folks have (of course!) done a ton of optimisation on two key things:

The DNA cleanup
The library preps

I think that for most DNA extractions cleanups can vary a lot, but something like a Blue Pippin will do an exceptional job if you have access to one and can use it (they won't work for SUPER long DNA fragments which can't move through a gel of course). The Blue Pippin also does your size selection, so you don't burn out pores sequencing lots of little stuff. This should help yield too.
Then it's all down to the library prep, and @gringer's answer has some good tips on that.

I am working with a fungus and am using CTAB method for DNA extraction and I used rapid kit for my initial library prep but it did not work and since than I have been using ligation kit and it has been giving me very good results. I managed to get data ranging from 4GB to 13GB.


text:
How can I systematically detect unknown barcode/adapter sequences within a set of samples?
I have often downloaded datasets from the SRA where the authors failed to mention which adapters were trimmed during the processing.
Local alignments tend to overcome this obstacle, but it feels a bit barbaric.
fastQC works occasionally to pick them up, but sometimes fails to find the actual adapter sequences.
Usually, I ended up looking up the kits they used and trying to grep for all the possible barcodes.
Is there a more robust/efficient way to do this?

This doesn't answer your question, but I hope there is a possibility to report such problems to SRA so that they ask the authors to publish the missing information.

Why do you feel that local alignment is a bit barbaric? It should be the default method in this day and age, unless you’re working with smallRNA sequencing. I tend to trim adapters to be on the safe side but I did lots of work without bothering, and just relying on local alignment.

The minion utility from the kraken/reaper toolkit may be helpful for this: http://wwwdev.ebi.ac.uk/enright-dev/kraken/reaper/src/reaper-latest/doc/minion.html

This looks like exactly the right type of tool. Although too bad it was designed mainly for the 3' end adapter. I wonder if you could just flip all your reads and apply it to the 5' end.

You mention that FastQC "fails to find the actual adapter sequences" - I guess you mean in the Adapter Sequence Contamination plot. However, the kmer and Sequence Content Plots are often useful even when the former fails. I've used these in the past - you can sometimes just read off the adapter sequence from the start of the Sequence Content Plot (or at least see how many bases to trim).

I'm not aware of any existing methods to do this, but here are a couple of ideas about how it might be done:
Canu has a method of adapter trimming which involves looking for the absence of overlap for reads. If there are no other reads which share sequence across a particular region, then the read is broken up at the point of low coverage, and small pieces are discarded. It would be possible to use a method like this to hunt for possible adapter/barcode sequences by preserving the short reads.
Another option is to do a kmer search at the start of reads, and see if any of the high-abundance kmers can be assembled together and/or matched to existing known adapters or barcodes.

If you happen to know a sequence that should be highly abundant in the library, you can grep its beginning or end (with pattern match highlighting) and see if the same sequence systematically comes just before or just after respectively. This kind of visual inspection can help you finding the adaptor.
For instance, in a previous lab, we were working on D. melanogaster small RNA sequencing data and my colleague knew from previous experience with this kind of data that the following small RNA was likely to be abundant: http://flybase.org/reports/FBgn0065042.html
We just had to grep it in the fastq file to see many lines with this sequence, next to another sequence that happened to be always the same: the unknown adapter.

May I know the reason of the downvote?
I've seen this method applied in a case of small RNA-seq, where one highly abundant sequence was expected. Visually inspecting the output of grep of this sequence (with pattern highlight) gave a very good hint of what the adaptor was (the non-highlighted part).

The question is asking about how to detect unknown adapter sequences, so the OP will not know about abundant sequences in advance. That's kind of the point of the question...

@tallphil I don't see the link between not knowing the adapter and not knowing of an abundant sequence expected to be present in the data. If I remember well, in the example I mention in my comment, my colleague knew from previous experience with this kind of data that the following small RNA was likely to be abundant: http://flybase.org/reports/FBgn0065042.html
We just had to grep it in the fastq file to see many lines with this sequence, next to another sequence that happened to be always the same: the unknown adapter.

Actually I just re-read your post and now I see what you meant. This is a reasonable idea. However, I think you explained it poorly in the sense that a reader might be confused and think you meant that searching for the most abundant sequence might come up with the barcode. You should have specified that the "abundant sequence" in this case was a known nucleic acid sequence that would be expected to have adapters ligated to one or both ends.

Ah yes, apologies - this is exactly how I read it. Probably didn't help that the original question mentions grepping for expected adapter sequences, so this was fresh in my mind :) Sorry @bli! Downvote wasn't from me so I can't retract it I'm afraid.

I tried to clarify my explanations.


text:
Which tools can detect chimeric RNA (fusion genes) from WGS or RNA-Seq data?
Given WGS data or RNA-seq data, which tools can I use to detect gene fusions?

Most of these use RNA-seq data, some use WGS data, and some use both. They are listed alphabetically. I will add to the list when I discover more.

Barnacle: http://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-14-550
Bellerophontes: http://bioinformatics.oxfordjournals.org/content/28/16/2114.long
BreakDancer: http://www.nature.com/nmeth/journal/v6/n9/abs/nmeth.1363.html
BreakFusion: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3389765/
BreakPointer: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3561864/
ChimeraScan: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3187648/
Comrad: http://bioinformatics.oxfordjournals.org/content/27/11/1481.long
CRAC: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4053775/
deFuse: http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1001138
Dissect: http://bioinformatics.oxfordjournals.org/content/28/12/i179.abstract
EBARDenovo: http://bioinformatics.oxfordjournals.org/content/early/2013/03/01/bioinformatics.btt092
EricScript: http://bioinformatics.oxfordjournals.org/content/28/24/3232
FusionAnalyser: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3439881/
FusionCatcher: http://biorxiv.org/content/early/2014/11/19/011650.full-text.pdf+html
FusionFinder: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3384600/
FusionHunter: http://bioinformatics.oxfordjournals.org/content/27/12/1708.long
FusionMap: http://bioinformatics.oxfordjournals.org/content/27/14/1922
FusionQ: http://www.biomedcentral.com/1471-2105/14/193
FusionSeq: http://www.genomebiology.com/2010/11/10/R104
IDP-fusion: http://nar.oxfordjournals.org/content/early/2015/06/03/nar.gkv562.full
iFUSE: http://bioinformatics.oxfordjournals.org/content/29/13/1700.long
InFusion: https://bitbucket.org/kokonech/infusion/wiki/Home
INTEGRATE: http://www.ncbi.nlm.nih.gov/pubmed/26556708
JAFFA: http://www.genomemedicine.com/content/7/1/43
LifeScope: http://www.thermofisher.com/no/en/home/life-science/sequencing/next-generation-sequencing/solid-next-generation-sequencing/solid-next-generation-sequencing-data-analysis-solutions/lifescope-data-analysis-solid-next-generation-sequencing/lifescope-genomic-analysis-software-solid-next-generation-sequencing.html
MapSplice: http://www.ncbi.nlm.nih.gov/pubmed/20802226
MOJO https://github.com/cband/MOJO
nFuse: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3483554/
Pegasus: http://www.biomedcentral.com/1752-0509/8/97
PRADA: http://www.ncbi.nlm.nih.gov/pubmed/24695405
ShortFuse: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3072550/
SnowShoes-FTD: http://nar.oxfordjournals.org/content/39/15/e100
SOAPFuse: http://www.genomebiology.com/2013/14/2/R12
SOAPFusion: http://www.ncbi.nlm.nih.gov/pubmed/24123671
STAR: http://bioinformatics.oxfordjournals.org/content/29/1/15
STAR-Fusion: https://github.com/STAR-Fusion/STAR-Fusion/wiki
TopHat-Fusion: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3245612/
TRUP: http://www.genomebiology.com/2015/16/1/7
ViralFusionSeq: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3582262/

Other useful programs:
Chimeraviz (visualization tools for gene fusions): https://www.ncbi.nlm.nih.gov/pubmed/28525538 (disclaimer: I created this)
Chimera: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4253834/
OncoFuse: http://bioinformatics.oxfordjournals.org/content/29/20/2539.long
FuMa: http://bioinformatics.oxfordjournals.org/content/early/2015/12/09/bioinformatics.btv721.abstract

1- Chimerascan + Star
2- Tophat Fusion
Two pipelines I have used extensively and recommend. The others listed by @L42 I can't speak to but Star Fusion sounds promising. Star is fast.


text:
How to transfer gff annotations in genome with extensive duplications?
Microbial genomes can contain extensive duplications. Often we'd like to transfer annotations from an annotated species to one that is newly sequenced.
Existing tools (e.g. RATT, LiftOver, Kraken) either make specific assumptions about how closely related the species are or fail to transfer when multiple matches are found in the new genome, especially if the sequences are highly similar.
Specifically, I have a synthetic biology application where genes can duplicate extensively. They are identical in sequence but duplicated many times and be relocated (i.e., not just adjacent to each other). None of the above mentioned tools are able to transfer coordinates of annotations to genomes with multiple copies of features.
Are there any pre-existing tools or software that transfer annotations in this scenario? Ideas for ways to do this robustly?

Can you provide more detail? What's your input (e.g. raw reads, assemblies, open-reading frames)? what kind of sequencing? what is the exact output you want? ' I'm not sure I understand your point about the assumptions. Also, LiftOver and Kraken are completely different tools with different uses

Would this even be theoretically possible? How can you assume the annotations are transferable if there are extensive duplications? It's probably better to look for homologs instead.

@terdon do you mean orthologs? homologs = orthologs (non-duplicated) + paralogs (duplicated)

@Chris_Rands No, I meant homologs. Precisely because we can't know whether they are ortho- or para- (I have a nice post on the difference between the two here, by the way), so all you can do first is find homologs and then try and figure out whether they are similar enough to carry any annotations over.

@terdon I see well resolving orthologs/paralogs is not easy of course, but it can be done, depending on the exact data (I don't know what the OP's data looks like), for example some of my colleagues maintain orthodb http://www.orthodb.org/

Oh, of course it can be done! My point was that looking for regions of homology (of whatever type) seems like a better way of transferring annotations than attempting to translate genomic coordinates between genomes of different species.

@Chris_Rands: Input would be assemblies, e.g. de novo from gDNA sequencing. Output would be a transfer of annotations (e.g. gff format) from characterized species to newly assembled genome (coordinate transfer). Both LiftOver and Kraken (this one, just to make sure were on same page: https://github.com/nedaz/kraken) do this. LiftOver more appropriate for coordinate transfer between closely related sequences e.g. different assemblies; Kraken uses genome alignment (MUMer, Satsuma) so better for more divergent sequences.

@terdon: distinguishing between type/origin of homology would go beyond scope of what I would want to accomplish, but the difference is important as you point out. also right to say that transferring smaller homologous regions would be better, especially for diverged species. fyi about application: I have a syn bio application where genes can duplicate extensively. they are identical in sequence but duplicated many times and relocated (i.e. not adjacent). None of the above mentioned tools was able to transfer coordinates of annotations to genomes with multiple copies of annotation.

Yes, nor would I expect them to. That's what I was saying. The liftover tools simply map coordinates, they won't be able to deal with this sort of thing. I am afraid you will have to do it manually by getting a list of genes/proteins of interest, finding their homologs and transferring the annotations over (with the obvious caveats about whether or not you can assume the annotations are transferable). Won't be much fun, unfortunately.

apologies I thought you meant kraken: http://ccb.jhu.edu/software/kraken/, who names these tools? anyway, this is quite non-trivial to do properly. you'll need to do genome assembly, gene predictions and ortholog/paralog assignment; there are various pipelines (some reviewed here: https://www.ncbi.nlm.nih.gov/pubmed/27043882), but they'll take some time. alternatively, for something more 'quick and dirty', @terdon 's suggestions seems sensible

There is one very simplistic way I use which might work for what you are doing, it is similar to what terdon proposed.
Take a de-novo microbial genome annotation tool (I have my own, but you could use/modify prokka). Tools like these often first predict gene boundaries (with other tools like prodigal or glimmer) and then try to assign a function to found genes. This function assignment is often done with BLAST and other tools ... and that is where you can go in and modify to do what you need.
I use a "knowledge" protein database of genes I want to have very strictly annotated as a first line of annotation (e.g. in your case: the annotated genomes). For that I loop through very strict identity/similarity parameters which get gradually relaxed.
E.g.:
Loop 0: only transfer annotations at 100% DNA identity, same length.
Loop 1: only transfer annotations at 100% similarity, same length.
Loop 2: only transfer annotations at 99% similarity, length +/- 1%.
...
Loop n: only transfer annotations at 100-(n-1)% similarity, length +/- (n-1)%.
In each loop, obviously only annotate what has not been annotated in previous loops.
After that, use "normal" annotation pipeline of the tool to annotate the rest.

Doesn't that require the target genome's genes to have been found first? Or  can your tool also do de-novo gene prediction? (sounds like a very useful tool, by the way, kudos!)

Prokaryotik gene finding/prediction is a more or less solved problem, existing tools work reasonably well. See http://prodigal.ornl.gov/ and http://prodigal.ornl.gov/ (just to name two).

Yes, I know, I was just surprised you didn't mention that in your answer. If I understand correctly, the first step would be for the OP to find the list of putative genes in their newly sequenced genome, right?

Correct. Prokka (http://www.vicbioinformatics.com/software.prokka.shtml) uses a whole battery of third-party tools (including prodigal) to annotate a genome de-novo, I started by modifying prokka before I wrote my own one, which uses some ideas from the prokka pipeline.

I think you will have to first identify the regions homologous to the ones defined in your GFF and then transfer the annotations. Of course, the assumption there is that the homolog will also have the same annotation which is often not true. However, I don't see how you can do it in any other way since you cannot use genomic coordinates (and you would still be making the same assumption even if you could, anyway) when the genomes are so different.
For a very simplistic approach (which might be enough if, as you say, your sequences are almost identical), you can do something like:

Collect the sequences of interest from your already annotated species.

Use a tool like genewise or exonerate to map these into the target genome. Both tools can return gff-formatted output and both can find multiple hits in the target genome. For what you want, I would suggest using a very high threshold of sequence similarity and query coverage (where the target sequence found covers all or most of the query sequence used).
Since these are microbial genomes and therefore splicing isn't a problem, you could do the same thing with even a simple BLASTn or tBLASTn if you start from protein sequences.

At this point, you should have a list of homologs (some of which will be orthologs and others paralogs) and you can transfer the annotations of the query sequence over to the target.

Again, I stress that this is making a whopping huge assumption: that homologous sequences have the same function and can automatically be annotated as whatever you had in the query genome. This is going to be true for many cases but it will also be false for others. Especially if you are looking at paralogs (genes whose duplication occurred after the speciation event and are therefore likely to have diverged in function).
However, as I said before, this problem would be exactly the same even if you did manage to transfer annotations just by identifying the syntenic regions of the genomes1, so there's not much difference there.

1 As I said in the comments, I don't see how this could be possible. By definition, if you have extensive duplications, the genomic coordinates will be completely different and it is impossible to map from one genome into the other.


text:
Rapid metagenomics classifiers on long read data
I recently used the minION (Nanopore, 9.4 flow cell, RAD001 kit) to generate a metagenome out of environmental samples.
Passed reads weren't brilliant (196, average 1,594bp lenght), but working with centrifuge the classification outputs turned out to have quite low hitLength to queryLength scores (average 2%, max 14%). This plus score values don't give me a lot of confidence towards the results I got.
Has anyone else used centrifuge and experienced the same?

Sorry, but "has anyone else seen this" isn't on topic here (that's just a discussion and asking for personal opinion). Could you [edit] your question to make it request a solution to a specific problem instead? Something like "what sort of threshold should I use for centrifuge" or similar?

Yes, this question can be improved, but I think people can understand what OP really wants. I don't think it deserves a closing.

You should also keep in mind that single read accuracy of ONT still is a bit lacking, a 2D or 1D^2 accuracy of 95% still means that on average there's one error every 20 bp, and due to the nature of the data some stretches may be more junky than others. Maybe centrifuge doesn't like that.
Getting a feel of what the data can tell you (and what not) is pretty important.
You may want to try other classification tools to get an idea how well centrifuge works, Kraken comes to mind. But also: what about BLASTing some of these reads at NCBI and simply look at the results to get a first impression?


text:
Changing the record id in a FASTA file using BioPython
I have the following FASTA file, original.fasta:
>foo
GCTCACACATAGTTGATGCAGATGTTGAATTCACTATGAGGTGGGAGGATGTAGGGCCA

I need to change the record id from foo to bar, so I wrote the following code:
from Bio import SeqIO

original_file = r"path\to\original.fasta"
corrected_file = r"path\to\corrected.fasta"

with open(original_file) as original, open(corrected_file, 'w') as corrected:
    records = SeqIO.parse(original_file, 'fasta')
    for record in records:
        print record.id             # prints 'foo'
        if record.id == 'foo':
            record.id = 'bar'
        print record.id             # prints 'bar' as expected
        SeqIO.write(record, corrected, 'fasta')

We printed the record id before and after the change, and get the expected result. We can even doublecheck by reading in the corrected file again with BioPython and printing out the record id:
with open(corrected_file) as corrected:
    for record in SeqIO.parse(corrected, 'fasta'):
        print record.id                  # prints 'bar', as expected

However, if we open the corrected file in a text editor, we see that the record id is not bar but  bar foo:
>bar foo
GCTCACACATAGTTGATGCAGATGTTGAATTCACTATGAGGTGGGAGGATGTAGGGCCA

We can confirm that this is what is written to the file if we read the file using plain Python:
with open(corrected_file) as corrected:
    print corrected.readlines()[0][1:] # prints 'bar foo'

Is this a bug in BioPython? And if not, what did I do wrong and how do I change the record id in a FASTA file using BioPython?

Keep in mind that record.id is not equal to the full record header, as record.description is a second part of it. If you parse the headers and they contain two space separated entries, record.id will only print the first one.

In addition, it seems you open the original file unnecessarily since you never do anything with the original object.

I think @Jenez is completely correct; Also you asked the same question on Biostars and got the same answer https://www.biostars.org/p/156428/

@Jenez why not post that as an answer?

@Chris_Rands yes, but we need content here so BioGeek was nice enough to ask again so we can have another answered post! Now the next person who comes here with a similar issue will find an answer. Yay!

Don’t write paths as r"foo\bar", this only works on Windows. Write them as "foo/bar" (no raw string needed), this works everywhere (including Windows).

@terdon Because I felt it didn't accurately answer the question fully. It was more of a nudge in the right direction.

Regarding file paths, I suppose it is even better to use os.path.join, which automatically choses the correct separator.

If I use SeqIO.parse(filehandle, 'fasta') to parse a FASTA file, then it will return a SeqRecord object where the id and name are the first word (everything before the first whitespace) of the line beginning with > and the description is the complete line (all not including the initial >). (This behaviour can overruled by providing a custom title2ids function).
So, if I for example have a original.fasta file like:
>Lorem|ipsum>dolor sit amet
GCTCACACATAGTTGATGCAGATGTTGAATTCACTATGAGGTGGGAGGATGTAGGGCCA

Then I will get:
>>> from Bio import SeqIO
>>> path = r'C:\path\to\original.fasta'
>>> records = SeqIO.parse(open(path), 'fasta')
>>> record = next(records)
>>> record.id
'Lorem|ipsum>dolor'
>>> record.name
'Lorem|ipsum>dolor'
>>> record.description
'Lorem|ipsum>dolor sit amet'

So, what happens in my case when you do
    if record.id == 'foo':
        record.id = 'bar'

is that the record.id is successfully changed from foo to bar, but that the record.description is not changed and stays foo. That's why, when the FASTA file is printed out, I see the described behaviour of
>bar foo
GCTCACACATAGTTGATGCAGATGTTGAATTCACTATGAGGTGGGAGGATGTAGGGCCA

So, the solution to my problem is to both change the id AND the description:
from Bio import SeqIO

original_file = r"path\to\original.fasta"
corrected_file = r"path\to\corrected.fasta"

with open(original_file) as original, open(corrected_file, 'w') as corrected:
    records = SeqIO.parse(original_file, 'fasta')
    for record in records:
        if record.id == 'foo':
            record.id = 'bar'
            record.description = 'bar'
        SeqIO.write(record, corrected, 'fasta')

with open(corrected_file) as corrected:
    records = SeqIO.parse(corrected, 'fasta')
    for record in records:
        print record.id             # prints bar
        print record.name           # prints bar
        print record.description    # prints bar


text:
Why is SMALT better for microbial genomics than other mappers?
SMALT seems to be one of the most used read mappers for bacterial data, see, e.g., this query. I do not say that it is not a great mapper, but I cannot easily see what are its main strengths compared to mappers such as BWA-MEM, Bowtie2, NovoAlign or GEM. Moreover, it is not even published.
Could you name some of its distinguishing features (e.g., user support by Sanger Pathogens)?
So far I have heard only arguments like "We use SMALT because everyone does it.", but this is not convincing enough for me.

Is is really that popular? e.g. compare your query to this one: https://scholar.google.com/scholar?as_ylo=2016&q="Bowtie2"+bacteria I suggest you re-phrase to just ask about the differences between SMALT and the other algorithms

There are big clusters of researchers who do not use anything else than SMALT. For instance, a lot of people in bacterial population genomics.

Yes, or this query. In fact, if we go by the results of the three queries, it looks like SMALT is an order of magnitude less common than either Bowtie or BWA.

@terdon BWA-MEM has 617 results. SMALT 162. I agree that it is less but it is still a highly used mapper and in certain communities the dominant one. Your last query (BWA) has a lot of results like biological warfare agents, etc.

Sure, counting google results is a horrible way of doing this. My point is that I don't see much evidence that this aligner is all that popular. I hadn't even heard of it, for example, and I've been working on NGS for the past couple of years and around people working on NGS for more than 15. I think you just happen to be in a community where this tool is over represented. I would have told you the de facto standard is BWA instead. I guess it just depends on what one has been exposed to.

I completely agree – I am now in a community where SMALT is highly overrepresented. I also consider BWA-MEM to be a standard (it is a little bit unfair for other mappers to talk about BWA as a whole since it comprises three distinct mappers). This is why I am trying find some reasons why SMALT is so popular in certain communities and to understand its distinguishing features from "more standard" mappers.

Heh, I refer to BWA because where I work we've found that bwa sampe + bwa samse gives better results than bwa mem :)

@user172818 Thank you!! This is exactly what I wanted to know.

Btw. What is interesting about the up/down-votes here is that whenever this question got 1 like, someone immediately put 1 dislike so the sum was 0 most of the time.

@user172818 could you please put your flexibility explanation in as an answer to this question; it seems more appropriate as an answer rather than a suggestion for improvements to the question.

It is also recommended mapper for reapr, because it is able to align mate pairs independently to each other (-x option). I do not really know if other mappers do provide this option as well, I have not investigate that...

I have used smalt a couple of times. Smalt is more flexible than bwa/bowtie and can be tuned to be more sensitive to divergent hits, which is useful to certain applications. I heard Sanger did evaluate a few mappers for some non-typical applications (not sure what). They found smalt to perform better. Also smalt was developed at Sanger. I guess they could modify smalt based on their needs, though I don't know the details.

Thank you for your answer. Btw, Do you known how SMALT compares to NovoAlign? My guess is that Novoalign could perform better, but I have no evidence for it.


text:
Gene set enrichment analysis on differential phosphorylation sites
I have:

A list of differentially phosphorylated sites in a knockout condition. Some genes contain as many as 70 possible phosphorylation sites; others contain only one.
A list of genes belonging to a specific gene set annotation.

How can I test the differentially phosphorylated proteins for enrichment of this annotation?
A few ideas I’ve considered:

Ignore the number of phosphorylation events detected within a gene and simply count the gene as differentially phosphorylated if it contains at least a single site that is differentially phosphorylated. Compare the enrichment score for this selected set to the enrichment score for the set genes containing at least one site that is not differentially phosphorylated. The problem here is that genes with a large number of phosphorylation sites have almost no influence in the enrichment score, since it’s almost certain that they’ll have at least one site that is differentially phosphorylated and at least one that is not.
Mark each candidate phosphorylation site as either “in the set” or “not in the set” based on the protein in which it’s found. Then perform the enrichment analysis using the set of annotated phosphorylation sites instead of using the traditional enrichment analysis performed at the gene level. The potential problem with this approach is that it may place too much influence on genes with many potential phosphorylation sites.
Aggregate all candidate phosphorylation sites within a gene and use some numerical threshold to determine whether the gene is differentially phosphorylated or not. (There are various ways that this could be done.) Then perform enrichment analysis using the resulting set of differentially phosphorylated genes. A possible problem here is that some of the phosphorylation sites may be more functionally important than others, so it’s not clear how to weight the relative importance of individual phosphorylation sites within a gene.

I realize the goal here isn’t well-defined mathematically; I’m mainly curious what approach makes the most sense given the biological context. Of the above approaches, I’m currently leaning towards approach 2 because it’s straightforward to implement and it at least attempts to account for the variable number of phosphorylation events within a gene.
NOTE: I also have normalized phosphoprotein and protein abundances for all of these sites, obtained from mass spec. So if the solution requires an alternative method of computing differential phosphorylation, that’s fine.

Your dataset looks perfect for the SUMSTAT enrichment test.

You need to come up with a statistic representing your gene. The simplest ideas here are the number of sites or, probably better, a proportion of phosphorylated sites.
Now you can compute a statistic for every gene set, for example just sum (SUMSTAT). You can also have average, sum of squares or something else.
Geta a null distribution of your gene set statistics by permutations. You can either randomize statistics for genes, or just randomly assign phosphorylated sites across the genome keeping the total number of sites.

Now you can compute p-value just by comparing your value with the null-distribution.
You have to be cautious of two things:

You are testing multiple gene sets, i.e. you perform a statistical test for every gene set. Therefore correct your significance for multiple testing. I suggest using FDR for this.
Beware of potential biases when doing permutations. The obvious one is gene length, i.e. longer genes have higher chance of getting at least one site. But there can be other ones like GC-content, chromosomes, etc. You can overcome this by using more realistic permutations, or by controlling for the potential correlations. You can get a few ideas on fighting biases from this paper (sorry for the self-advertisement).


text:
How do I validate a single sample ArrayCGH result?
We have arrayCGH (aCGH) results for one sample. There is a 0.5 Mb terminal duplication on chromosome 19 (62995490-63407936, according to NCBI36/hg18). The duplication is rare: a literature review suggests there are only 3-4 samples with clinical information.
What are the steps to validate the results? How do we ascertain that this duplication is the cause of the clinical symptoms?
I have some ideas:

aCGH the parents. Not sure how this would help.
whole genome exome sequencing. Worried this might make it more difficult to pinpoint genetic cause.
whole genome sequencing?
other ideas?

Note: I am new to aCGH and high-throughput sequencing, any advice is welcome.

Is the duplication in the telomere or before it?

@DevonRyan I think, it is before, on the last band of chr19. Does it matter, sorry not being ignorant, just no clue about it.

Telomeric sequences are highly repetitive, so yes it matters :)

Telomere starts at the end of the last band, or is it part of chromosome?

If you look at the UCSC genome browser, the telomere is the GGGTTA simple repeat at the end. In others it's just a stretch of NNNN in the reference (no clue how this looks in hg18).

Just as a little suggestion, especially in the beta you should take the recommendations of the SE interface very seriously and always leave a comment when down-voting, or you demonstrate just the prejudice of hostility! Also be a little open minded, questions about omics data and how to validate omics results are typically very relevant. Otherwise you will notice that nothing can be asked anymore. Just saying, from and old bs mod. P.S.: I have the feeling that people down-vote because they do not understand, that's not always the fault of the Op.

Another observation, this question might be ethically questionable! You are putting forward info on a real case to diagnosis, possibly identifying information is revealed.

Why not just good old qPCR? That's (A) quick, (B) cheap and (C) easy to analyze. If you care about the exact location of the break point (I'm guessing from the context that you don't), then targeted sequencing with a custom capture kit would work.
Regarding validating the biological relevance of this, there are multiple (parallel) routes one can take. Firstly, screen unaffected family members for this same alteration. If you find this in unaffected individuals then it's obviously not the causative alteration. Ideally, one would also do either a cell-line experiment or a mouse (or other model organism) experiment to see if you can reconstitute at least some component of the clinical phenotype. This may not always be possible, of course.


text:
Filtering imputed GWAS SNPs based on a MAF difference of 10%
There are many posts on the web regarding QC steps pre and post-imputation.
Does applying below (new?) 10% MAF difference rule make sense, pitfalls?
Here is the process:

Get MAF for imputed set, using SNPTEST with flag -summary_stats_only
Convert imputed set to hard-calls using gtools with flag --threshold 0.9
If the MAF from step 1 and step 2 differs more than 10% than exclude the variant.

More info:

GWAS is 50K vs 50K case control samples.
This step is applied after info > 0.4 filter.

I cannot think of any principled rationale for choosing this filtering strategy.
However, I am going to take a guess that this filtering strategy is supposed to filter out SNPs for which imputation did not work well?  In that case the appropriate statistic to filter on is the INFO score as described here.
You might consider picking a higher threshold than 0.4 if you feel that you aren't filtering out enough variants.


text:
Working with old genome builds
Is working with and relying on old genome builds still valid?
For example NCBI36/hg18. Would results from papers based on old builds require LiftOver and re-analysis to be useful?
A bit of context, this is related to other post, where we have aCGH results based on old build: How do I validate a single sample ArrayCGH result?

This will be probably dependent what type of analysis you have in your mind. In the end, all the data we generate today will be one day obsolete, but it does not necessary mean that all the conclusions are wrong. If you would be more specific about types of analysis you have in your mind (or concrete papers using hg18), maybe it would be easier to give a correct answer.

In my opinion, it is not very reliable. LiftOver is very limited in terms of transformations it can support. The LiftOver Chain format can capture only matching regions in the same order. It means that it can account for indels, but even simple structural variations become problematic.
For instance, when a newer assembly is available, it is usually a recommended practice to remap all the reads rather than transform the existing alignments.

I think that right now, the only human builds that are worth considering are hg19/GRCh37 as many data bases such as gnomAD still exclusively use this release. On the other hand, hg38/GRCh8 has many important fixes and the useful (but yet underused) feature of alternative loci.
Anything from older releases should be remapped to a more recent one.

You could use liftOver which isn't always great.
Whenever I encounter this (especially NGS data readily available on the SRA), I often just get the raw files (e.g. fastqs) and re-align/re-map.
In your case (arrays) it may be a bit tough. Not impossible though, as I recently took some old yeast DNA/RNA microarray data and updated it to the newest genome. Just requires the right data (like DNA for normalization) and a good understanding of the entire process.
A last resort/alternative is to align your new data to the old genome to be able to make comparisons. This isn't ideal but works in cases where upgrading one source isn't possible or is a HUGE amount of time/effort. I've done this for a few fly experiments where all the available/previous data was done in dm3. All the old genomes can generally be found on http://archive.ensembl.org .

For mouse, I still see people using mm9/NCBI37 in high profile publications even though mm10/GRCm38 was released more than 5 years ago (2011). I personally don't think that's a great idea, but it's certainly valid according to the peer reviewers.
It also depends on your application. If you are working with coding regions (likely well known for a long time) or extracting genome-wide stats (enrichment at TSS, for example), the differences should be negligible.

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


text:
Is there a point in recalibration of scores for variant calling?
The most variant calling pipeline GATK include a Base Quality Score Recalibration (BQSR) which requires a list of known variants. Recently, some work has been done for reference-free recalibration of scores as well: Lacer and atlas, which is motivated by making the most for aDNA and low coverage datasets.
The importance for aDNA is explained in this lecture, but it is not clear to me if / how is important BQSR is for fresh DNA samples with decent (>15x) coverage. Especially when I work with non-model organisms and I can not simply use the standard tools.
How big an impact does recalibration of scores have on variant calling? Is there a rule of thumb for which it is / it is not worth the effort?

I personally don't think BQSR has a huge impact on variant calling, but you don't really need to guess. If you run GATK BQSR, it outputs a table and charts of exactly how much quality scores are adjusted. The adjustment will vary depending on the position in the read and genomic context (previous and following base). In my experience, the difference is a few points at most, but it's certainly noticeable.
GATK recommends BQSR for both genome and exome data, which is normally much higher than 15x.

That's a good question.
I'd say that you don't need to bother with variant recalibration for

low number of samples (e.g., just two trios); I could not get GTAK recalibration of variant scores to work anyway
high-coverage samples (e.g., X Ten genomes with 30x coverage) where the DNA samples themselves are of high, comparable quality and have been sequenced with consistent technology.

Generally, it is my impression that a lot of the thoughts and advanced statistical models built into GATK come from the earlier phases of the 1000 Genomes project. This means (1) low-coverage, (2) different coverage genomes (3) sequenced with varying technology versions by (4) different samples and (5) population sequencing.
If you are in a clinical setting where you do 30x sequencing on X Ten platforms only anyway, then variant recalibration will probably not help you that much.
On the other hand, if you are integrating many data sets from different data centers and machine versions etc., variant recalibration might be worth a shot.
A good check would be looking at genotype quality distributions and other variant/quality related metrics before and after recalibration.
Anyone: please correct me if I'm wrong!

Are you talking about base quality score recalibration (BQSR) here or about  variant quality score recalibration (VQSR)? I think the OP is referring to BQSR but you are discussing VQSR.

Yes, OP confirmed. The question is about BQSR so I am afraid you're answering the wrong question.

sigh and there I thought I could contribute something.

Ideally these BQSR methods were made keeping in mind how technical errors will actually screw up the base quality calls and when the machines were still more on development phase while being used for the 1000G project. As of now machines are more powerful and strong where it will be unlikely to use it but still we use with listed SNPs to find the covariates and build a model around the data using the information with machine learning tricks to improve the quality of those base calls. Ideally it should be more appropriate when old machines from Illumina or other standard companies are being used but with new machines which are much powerful and having high throughput they should tend to go down. I do not recall if such tests have been made but obviously I know new sequencing machine always make such tests to show that they have reduced such errors but still recommend such BQSR for variant calls. Now the problem is the list of SNPs, this to me is the real problem since the list we use is far from being Gold standard and if that is not properly taken care of everything we infer about quality is still shaky. This link is pretty informative but it's an old one. I would really see improvements with new sequencers.  However very less people care about such tests in academic research and also translational lab will really not invest time and money on such unless the facility has some bioinformaticians who always does such testing while buying a new sequencer for the institute. In terms of clinical genomics for finding variants I reckon most powerful and up-to-date sequencers should be used but not sure if they still use BQSR and if so what is the list they use to build model of covariation around the data.

Note that the question is motivated by non-model organism research - I can not use a list of known variants, because such a list does not exist for my species. Therefore I wanted to know, how important is to recalibrate QS, because thechinally it is possible, but not just by running one more step in GATK pipeline.

I agree that it is non-model organism and that is why you will not have any such list of variants. But since the motivation of the approach was also asking about the BQSR so I said. You can take a look at this https://media.readthedocs.org/pdf/lts-workflows-sm-non-model-toolkit/latest/lts-workflows-sm-non-model-toolkit.pdf about how to use your HC variants of your sample to recalibrate. Also this gatlk link might help as well. http://gatkforums.broadinstitute.org/gatk/discussion/3286/quality-score-recalibration-for-non-model-organisms. Now the decision lies on your hand as to use or compare.

@KamilSJaron for non-model this is one way to do but if your sequencing machine is pretty new and improved accuracy you might also do away with the step. I would read publications to see what they do but still for my sake do calls without BQSR and with BQSR with HC variants and use them as database and compare to come to conclusion myself. Thats my opinion. It also depends on the virtue of the project as well.

The first link is very relevant to my original question: "While GATK UnifiedGenotyper suffers during indel calling without recalibration and realignment, both HaplotypeCaller and FreeBayes perform as good or better without these steps.". Thanks. The second link is relevant as well, but I do not have enough individuals sequenced to chose their approach for recalibration.

@KamilSJaron I am glad it is relevant but then again I would be saying that you can run both with and without and make some estimation. Since you have not much samples to create your own HC SNP database you can also do it with strict stringent SNPs from your individuals. Or since you do not have much samples just avoid the BQSR step and pull out top variants and not large fraction of variants. The top variants even though the score might not very accurate but the calls will be still high confident and preferably true positives. I guess depends on number of variants you stream down to.

BQSR is still relevant, one of the reasons for this is some new illumina hardware such as the NextSeq can only output binned Q scores, BQSR essentially "un-bins" the Q scores giving more granularity which has it's uses in deep somatic variant calling, http://gatkforums.broadinstitute.org/gatk/discussion/4594/beware-of-using-binned-quality-scores-with-some-gatk-procedures.  Additionally the NextSeq suffers from a high confidence poly G issue, BQSR will be beneficial here too https://sequencing.qcfail.com/articles/illumina-2-colour-chemistry-can-overcall-high-confidence-g-bases/

Ah now that was something I was unaware since I have not yet had to encounter with NextSeq. This is a good point. Then the OP should also think of it and yes I was merely considering how the BQSR in first came into picture. But this is a good catch for binned Q scores and about the poly G issue.

In case BQSR is not an option (i.e. non-model organisms) it would be best to use some internal control sequence such as PhiX for illumina platform. Although this is supposed to be common practice some facilities ignore it. In principle, the machines should use these sequences as reference so that the scoring would be more accurate. In my experience the first 10-15 bases of the illumina reads always had lower quality. This can easily be seen in the nucleotide distributions. I would advise trimming the first 10-15 bases and quality based end-trimmingiIf the quality of the individual reads are important, such as low coverage resequencing or de-novo genome assembly applications.


text:
Downloading a reference Genome for Bowtie2
How do I download a reference genome that I can use with bowtie2? Specifically HG19. On UCSC there are a lot of file options.

tl;dr: Just use the either the downloads on the Bowtie2 homepage or the Illumina iGenomes. Or just uncompress and concatenate the FASTA files found on UCSC goldenpath and then build the index.
A bit longer answer:
There are two components to "genome for a read mapper" such as Bowtie or BWA.
First, you need to choose the actual sequence (genome release such as GRCh37/hg19 or GRCh38/hg38). There are patch releases such as GRCh37.p3 where some bases might be exchanged and depending on the release, some "unmapped" loci contigs might be added, but generally GRCh37.p1 is roughly the same as GRCh37.p2, for example. Usually, people have agreed on some specific patch version for each read and use this for read mapping.
Generally, there is the UCSC flavour hg19/hg38 etc. and the NCBI/GRC flavour GRCh37, GRCh38 etc. (similar with mouse). UCSC has no versioning besides the genome release and (to the best of my knowledge) does not update the genome sequence after releasing a hg19 FASTA file.
Second, you have to build the index files for each genome. Depending on the read mapper you use, you might or might not need the original FASTA files for the alignment. For Bowtie and Bowtie 2, you don't need the original FASTA files after building the index as Bowtie 1/2 can reconstruct the sequence "on the fly" from the index files.
HTH

I don't know how I managed to miss the download on the bowtie homepage. Hope this helps someone else!

It’s a matter of preference I guess but I recommend the Ensembl builds. Decide whether you want the toplevel or primary assembly, and whether you want soft-masked, repeat-masked or unmasked files. The naming schema is very straightforward; the combinations are described in the README file, and all files reside in one directory.
For example, if you want the unmasked primary assembly, the file to download would be Homo_sapiens.GRCh37.75.dna.primary_assembly.fa.gz.
As for GoldenPath/UCSC, there’s no need to download and concatenate separate chromosomes (contrary to what the other answer said); you can download the whole (toplevel) reference from the bigZips directory; from the README:

This directory contains the Feb. 2009 assembly of the human genome (hg19,
  GRCh37 Genome Reference Consortium Human Reference 37 (GCA_000001405.1)),
  as well as repeat annotations and GenBank sequences.

There are essentially three options here:

chromFa.tar.gz, which contains the whole genome in one chromosome per file;
chromFaMasked.tar.gz, the same with repeats masked by N;
hg19.2bit, which is the whole genome in one file, but needs to be extracted using the utility program twoBitToFa, which needs to be downloaded separately.

In any case, I always download the reference and build my own index for mapping, since this allows me more control; not everybody might need this much control, but then building the index once is fairly fast anyway.

I think this triggers another question "what's the difference between different versions of the same genome build?". The question's answer should include the difference between DNA and RNA-seq/functional genomics analysis.
In the DNA/variant world, people will generally stick to whatever the large sequencing projects/Heng Li decides is "best". In the RNA-seq/functional genomics world, careful curation of genomes is important, depending on the read mapper and also what downstream tools support (larger set of tools means longer tail of less used tools having idiosyncratic requirements).


text:
Finding peaks and estimating cell population sizes in multi-dimensional flow cytometry data
Our research institute processes a lot of flow cytometry data, but the produced data is under-utilised due to the effort required to process it. A typical run will produce 5 million events (ideally one event per cell), with up to 14 dimensions of [ideally] log-normal fluorescence values for particular groups of cells ("populations"). Due to various systematic errors, negative values, scatter, and a non-zero "zero" value can happen, but I'm going to ignore those for the purpose of this question and assume that the data are well-distributed.
Researchers will typically probe these data using manual filters (set up specifically for each experiment) to find populations of interest. I suppose a picture might help. This one shows three identifiable cell populations, in order of size one at about (X/Ly6C:2,Y/CD86:2), one at about (4.5,3), and one small population at (2,4).

Here's another plot that has two cell populations that are close to each other, such that the population "humps" overlap substantially.

Manual filters are typically used because it can be very difficult to distinguish between a noisy data point from a large cell population and a less noisy data point from a small cell population, particularly when considering populations that make up about 0.01% of the total cells.
As an additional complication, counting these populations can be difficult when populations overlap (as in the second image). A filter/slice through the plot that separates populations could count many cells as being members of the wrong population.
Can these populations be detected in an automated way? If it is assumed that the populations are spread in a gaussian fashion at some point in each dimension, is there some method that can be used to approximate the number of cells in each population, even when populations are close by?

There are quite a few algorithms developed for the automatic classification of multidimensional flow cytometry data, you can see a (not so recent?) review here: http://www.nature.com/nmeth/journal/v10/n3/full/nmeth.2365.html. In your case, you are interested in the unsupervised methods, since you do not have data from single populations (which would provide you "labels" associated to each sample), as I understood.
In other words, if you have your data set as a matrix of 5 million samples (cells) with 14 dimensions (FC values) each, then you are looking basically for algorithms that will do a clustering of the samples into groups. By the way your description of each sample belonging to a multivariate Normal distribution leads to a Gaussian Mixture Model, which can be estimated from your data using e.g. the scikit-learn library for python: http://scikit-learn.org/stable/modules/mixture.html

Automated gating methods are gaining popularity over manual and time consuming analysis in say FlowJo, you should take a look at the openCyto bioconductor package this is a framework which builds on top of existing Bioconductor infrastructure for flow cytometry (PLOS compbio paper) with the aim of making various analysis routes more accessible to the researcher by assisting getting data into and out various gating algorithms. You can find a nice intro here which starts with importing a manual gating scheme from FlowJo and a how-to on 1D and 2D gating can also be found here.
An excellent YouTube intro to analysing flow data in R by Ryan Duggan Cytometry on Air is worth watching especially for those which have no R experience and are curious as to what can be achieved.  Additionally it's worth having a look at the growing list of Bioconductor packages and datasets that exists for flow data: https://www.bioconductor.org/help/search/index.html?q=cytometry/


text:
Detecting structural variants with MinION data
Working on various cancers I have an interest in detecting structural variation (SV) in human, we've successfully used various tools like Pindel, SVDetect, Manta, and LUMPY, to name a few for detecting SVs in illumina short-read sequencing.  I curious if anyone has successfully used ONTs MinION sequencer for detecting SV, as there are many cases where longer reads would be beneficial for SV detection especially where events are long and/or occur in repetitive regions.  Has anyone tried or had success with tools previously designed for PacBio data such as Sniffles?
Supplementary question, obviously the MinION has not quite got the throughput of say the PacBio SEQUEL so can't generate low-pass coverage for the whole human genome as easily, has anyone got experience trying to generate reads targeting specific areas of the genome in areas that are known for translations, long inversions, duplications, etc?

Since I last looked, it appears Sniffles now explicitly mentions MinION data as an input on it GitHub page.  So that appears to be one tool which has been tried.

What's your coverage?

We're planning to try a targeted approach in the first instance, but we're currently at planning stage, so I can't give you a figure as yet.

There was a Structural Variant breakout session at the London Calling conference this year. Unfortunately I didn't attend that session, but MinION community members have access to Constance Donnell's summary of that:
https://community.nanoporetech.com/posts/breakout-structural-varia
Here are my attempts at grabbing non-creative chunks from those notes:

Professor Wigard Kloosterman has been developing a bioinformatics pipeline called NanoSV for mapping genomic structural variants in patients with congenital abnormalities
Tomas Sesani at the University of Utah studied the precise copy number of duplicated genes in individual virus genes over time, and tracked the SNP within copy number variable regions of the genome [no public tool mentioned]
Dr. Sudha Rao from Genotypic Technology in India used the MinION sequencer with Sniffles for structural variation calling

ONT have said that the sessions/talks from London Calling 2017 will all be made publicly available at some point in the future.

Thanks for this, yes it's often hard to figure out what software is out there for MinION related things as it ends up being posted behind the closed doors that only MinION community members have access too.  This is an issue for bioinformatics people because the often won't have the account since they don't own the instrument.