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?

source: Stack Exchange

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

sort: votes


text:
Cicros plot issue using two data track
Im trying to plot log fold change for a condition as well as no of patient for a particular gene mutated information but whenever im trying to do that i see only only the fold changes labeling not the patient but my idea is to see the fold change as well as mutation information for that gene im posting the .conf file along with the data files im using any help would be highly appreciated .
Or there is a flaw in my idea to show both of them which im not sure but i guess whatever i seen so far examples it can work!
My conf file code and my data files data
use_rules = yes
karyotype = data/karyotype/karyotype.human.txt
# spacing = 10000000
#
# you can write
#
# spacing = 10u
#
# See ticks.conf for examples.

chromosomes_units = 1000000
#chromosomes_display_default = no
#chromosomes                 = /hs[1234]$/
#chromosomes=hs5;hs7; hs17; hs11; hs9

<<include ideogram.conf>>

<<include four_ticks_4chr.conf>>

<image>
<<include etc/image.conf>>
</image>
<<include etc/colors.conf>>

<<include etc/colors_fonts_patterns.conf>>

<<include etc/housekeeping.conf>>

<plots>
<plot>
#this is for two direction expression data visualisation
type = histogram
#file=data2/HSC_VS_CMP.txt
#file = /home/punit/circos-0.69-5/data_files/LEUKEMIA_value_tab.txt
file = /home/punit/circos-0.69-5/data_files/CIRCOS_VAL1.txt

r0   = 0.60r
r1   = 0.35r
thickness = 10
min=-1.5
max=1.5
extend_bin = no

<rules>
<rule>
condition    = var(value) > 1.5
color        = dred
fill_color   = dred
</rule>

<rule>
condition    = var(value) < -1.5
color        = dgreen
fill_color   = dgreen
</rule>
</rules>
</plot>
#########################################################

###########################################################

<plot>
type    = text
color = black
#file = /home/punit/circos-0.69-5/LEUKEMIA_gene_tab.txt
file = /home/punit/circos-0.69-5/CIRCOS_LOC1.txt

r0   = 0.75r
r1   = 0.95r
show_links     = yes
link_dims      = 0p,0p,50p,0p,10p
# link_dims      = 4p,4p,8p,4p,4p
link_thickness = 2p
link_color     = grey
label_snuggle             = yes
max_snuggle_distance      = 2r
label_size   = 25p
label_font   = condensed
padding  = 0p
rpadding = 0p
</plot>

###############################################################################

#############################################################################################

    <plot>
    #this is for two direction expression data visualisation
    type = histogram
    color = blue
    #file=data2/HSC_VS_CMP.txt
    #file = /home/punit/circos-0.69-5/data_files/LEUKEMIA_value_tab.txt
    file = /home/punit/circos-0.69-5/M_DATA.txt

    r0 = 2.0r
    r1 = 1.8r
    orientation = out

    #max       = 20
    #min       = -20
    stroke_type = outline
    thickness   = 18
    extend_bin  = yes
    fill_color = blue

    <rules>
    <rule>
    condition    = var(value) > 10
    color        = dblue
    fill_color   = dblue
    </rule>

    <rule>
    condition    = var(value) < 8
    color        = dyellow
    fill_color   = dyellow
    </rule>
    </rules>

    </plot>

    </plots>

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


text:
How to combine multiple kernels of large sample datasets?
I have multiple large sample datasets in matrix format (each has 15000 rows and 5-50 columns) corresponding to different experiments. Each matrix contains the same number of samples(rows) but the variables(columns) are not the same. My objective is to cluster the samples on the basis of all the experiments.
I tried to use Unsupervised multiple kernel learning (UMKL) to integrate the datasets followed by kernelPCA using "mixkernel" package in R (https://www.ncbi.nlm.nih.gov/pubmed/29077792). The UMKL step calculates kernels for each dataset and combines the kernels using 3 different approaches: 1) calculating a consensus kernel from multiple kernels 2) calculating a sparse kernel preserving the original topology of the data 3) calculating a full kernel preserving the original topology of the data
The kernel calculation step was fine but the kernel integration step (all 3 approaches) runs very long and my computer hangs.
Is there any way to handle this problem? More specifically, is there any way to handle multiple kernel integration for large sample datasets?
Any suggestion alternative to using kernel methods will also work.

Welcome to Bioinformatics Stackexchange @SD1024. A recently published algorithm known as MOFA (Multi-Omics Factor Analysis, paper, github) is generating a lot of interest, and is designed answer exactly the sort of biological question you are describing. It claims to extract axes from multiple matrices with overlapping samples, but not necessarily overlapping features and handle missing data elegantly. I havn't tried it myself you, but I'm quite excited to do so.

Thanks for the suggestion. Actually, I used that package. If I have understood correctly, MOFA extracts the latent variables from individual data modality nicely but such variables do not integrate multiple modalities together. For example, If multiple Chip-Seq and CpG methylation assays are done on multiple genomic sites and across multiple patient samples/cell lines, MOFA extracts one set of latent factors for ChIP-seq and one set of latent factors for CpG methylation. It does not extract any latent factor (or a set of latent factors) that represents both the assays simultaneously.

As I say, i've not actaully used it myself, but the paper definately claims that the factors intergate multiple assays. For example, see their figure 2B, where the plot shows contributions to factor 1 from all assay types, but only contributions from drug response, mutation data and transcription (i.e. not from methylation) to factor 2

You are right. It is described in equation 1 of materials and method. A great help. Thanks a lot. :)

You can use regularized general canonical correlation analysis to cluster those matrices. I would suggest RGCCA, which allows to do this for an arbitrary number of matrices but it  doesn't use kernels.
There has been some efforts to use kernels on the canonical correlation analysis but there isn't any implementation (yet), you can read more on the related literature, specially the Kernel Generalized Canonical Correlation Analysis.
However in this kind of analysis you will need to come up with a design of how those matrices are related between them. It is not like a PCA where you give some data and you get some data back that summarize the matrix.

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


text:
peakranger bcp segmentation fault
I'm running into a segmentation fault running peakranger's bcp on Arch Linux. Binaries were compiled from source just today. Any idea how I can prevent this?
$ peakranger bcp -d chip_seq/align_count/Sample_DNA_loA_h3k4me3/alignment/Sample_DNA_loA_h3k4me3.hisat2.bam -c chip_seq/align_count/Sample_DNA_loA_input/alignment/Sample_DNA_loA_input.hisat2.bam --format bam -o chip_seq/align_count/Sample_DNA_loA/call_peaks/bcp --gene_annot_file /data/sai/annots/musmus_GRCm38.94.gtf.gz --verbose
program version:          1.18
Data files:
 File format:             bam
 Sample file:             chip_seq/align_count/Sample_DNA_loA_h3k4me3/alignment/Sample_DNA_loA_h3k4me3.hisat2.bam
 Control file:            chip_seq/align_count/Sample_DNA_loA_input/alignment/Sample_DNA_loA_input.hisat2.bam
Qualities:
 P value cut off:         0.0001
 sliding window size:     500
 Read extension length:   200
Output:
 Regions:                 chip_seq/align_count/Sample_DNA_loA/call_peaks/macs2_region.bed
 HTML reports:            Enabled
 Plot region length:      6000
 Annotation file:         /data/sai/annots/musmus_GRCm38.94.gtf.gz

Reads statistics:

 Treatment reads +:       2125600
 Treatment reads -:       2113802
 Average read length:     20
 Control reads +:         7582218
 Control reads -:         7575995
 Average read length:     20
 Verifying reads...
Warning: JH584295.1 only contains reads in the positive strand. The chromosome is removed.
Warning: GL456382.1 only contains reads in the negative strand. The chromosome is removed.
Warning: No reads were found in GL456382.1 of the treatment dataset. The chromosome is removed.
Warning: No reads were found in JH584295.1 of the treatment dataset. The chromosome is removed.

Reads statistics after correction:

 Treatment reads +:       2125597
 Treatment reads -:       2113800
 Control reads +:         7582196
 Control reads -:         7575947

 Calling peaks...

Segmentation fault (core dumped)

Reads were trimmed with trim galore and aligned with hisat2.
Update 2018-10-29: additional information of input and machine
$ samtools flagstat chip_seq/align_count/Sample_DNA_loA_h3k4me3/alignment/Sample_DNA_loA_h3k4me3.hisat2.bam
33521569 + 0 in total (QC-passed reads + QC-failed reads)
946323 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
5185725 + 0 mapped (15.47% : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
$ samtools flagstat chip_seq/align_count/Sample_DNA_loA_input/alignment/Sample_DNA_loA_input.hisat2.bam
22508880 + 0 in total (QC-passed reads + QC-failed reads)
3359461 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
18517674 + 0 mapped (82.27% : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

I tried this with two other input bam files and also hit a segfault.
Machine is a 36-core (72 hyperthread) workstation with 500GB of RAM.
$ peakranger --version
1.18
$ g++ --version
g++ (GCC) 8.2.0
$ uname -r
4.18.3-arch1-1-ARCH
$ pacman -Q boost
boost 1.67.0-7
$ pacman -Q zlib
zlib 1:1.2.11-3

Thanks @Llopis, I've clarified the question and added information on the bam files and the machine.


text:
Compare alignment quality of multiple sequencing runs aligned against the same reference genome
I have run Oxford Nanopore Technologies' MinION sequencing on the same DNA sample using three flowcells, each aligned against the same reference genome (E.coli K12 MG1655) using both BWA MEM and GraphMap and stored as BAM files.
How can I quantitatively and efficiently analyse the quality of alignment (percentage identity, insertion rate, deletion rate) of each of these files?

Qualimap will do this for you.

Go to qualimap.bioinfo.cipf.es
Run qualimap (default params are fine) on each BAM file
Open up the HTML output, and you can read off the %identity (they measure the opposite, i.e. mismatch rate, but 100% - mismatch rate is %identity of course), indel rate, etc.

One thing to watch out for (you don't mention it in your question, but just in case) is that you cannot directly compare Q scores - these are a bit of a mess and calculated very differently in each piece of software.
Unsolicited suggestion: you might also try NGM-LR for mapping MinION data. We've found it beats the others for our data (though we map to a distant reference).

Thanks for the suggestion. To clarify, is this the aligner you were referring to? https://github.com/philres/ngmlr

Yep. Poster is here: http://schatzlab.cshl.edu/publications/posters/2015/2015.GI.NextGenMap-LR.pdf
You could check out Sniffles too - detects structural variants from 3rd gen data. https://github.com/fritzsedlazeck/Sniffles

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


text:
Visualisation of long read RNA-Seq splicing
I have a dataset of Oxford Nanopore cDNA reads. Many of my reads are full-length or close to full-length transcripts, and I and am interested in examining alternative splicing. For this, I would like to begin by visualising my reads and comparing variants qualitatively.
I have tried visualising my data in both IGV and SeqMonk, but neither has given a satisfying result.
IGV

IGV shows some links between aligned sections, shown as fine blue lines and thick black lines. I cannot work out what the difference between these lines is, and even more confusingly, we have numerous alignments where two aligned exons come from the same sequence but are not joined by any line. This means that, to confirm or deny a spliced variant, I need to manually examine the read ID of each exon.
SeqMonk

SeqMonk gives a much cleaner visualisation, but unfortunately shows no links between aligned exons, and worse, I cannot find the read IDs. This means there is functionally (as far as I can see) no way to tell which exons come from the same read.
EDIT: Sashimi Plot
As recommended in this answer, I have tried using IGV's Sashimi Plot. Unfortunately, I believe the high error rate and relatively low coverage is causing this to create somewhat confusing output. The Sashimi Plot shows mostly nearly every junction to have just a single read supporting it.

Is there a visualisation tool which allows simple examination of splicing of full-length transcripts?

Are you sure you are mapping the reads correctly? Do they have a lot of soft clips, perhaps suggesting that they aren't fully mapped?

I couldn't say for sure, but it looks okay to me.. They were mapped with GMAP, and with read lengths in the range of 1500-4000bp we have soft clipping sitting around 20-50 for most of those reads. About half of the reads have mapping quality 40 (the others are very low, about a quarter have MAPQ 0.)

IGB gives the intended result without much hassle. It doesn't have the bells and whistles of IGV, but presents a clean and intuitive view of the individual reads.
In IGB, each read is presented as a single, contiguous line.

If you zoom in you can see the read ID for each read.

For a more holistic view, IGV's Sashimi plot can be modified to exclude spurious junctions supported by few reads. Here is a section of the above gene with default settings:

The minimum junction coverage can be set by right-clicking on the Sashimi plot itself, or application-wide under View -> Preferences -> Alignments -> Splice Junction Track Options. With minimum junction coverage 3, the noisy plot above now looks like this:

The numbers are still hard to read due to numerous overlapping almost-identical junctions, but at least the plot is now readable.

DNASTAR's software is for purchase, but high quality. GenVision Pro does genomic visualization, including Sashimi plots.
Edit: not sure why this answer is being downvoted, unless it's because the software isn't free. OP has tried IGV and SeqMonk, I mentioned an alternative he might not have heard of.
Here is a video demonstrating the use of Sashimi plots in GenVision Pro:
http://www.dnastar.com/t-support-videos.aspx?video=YJvcERoSIsg

Please downvoters, explain the reason for downvoting.

Maybe add example plot outputs?

The Sashmi Plot feature built into IGV. It gives a nice summary of the spliced transcripts and the coverage of each exon.
For example:


text:
Split FASTQ and matching BAM into matching chunks
I am running a slow downstream analysis on a large set of nanopore reads (approx 3 million), and would like to split them into smaller chunks, run the analysis in massively parallel, and then recombine. Originally I just split the FASTQ into chunks, re-aligned each chunk, and then merged the output, but here I would like to use an existing alignment so I can compare results with existing analyses (i.e. the alignments must be the same).
How can I efficiently split a FASTQ file and a BAM file giving the alignment of the FASTA file into chunks, ensuring that all of the reads in FASTQ chunk 1 are in BAM chunk 1, vice versa and so on?
My FASTQ is approximately 45GB and my BAM is 33GB, so I would prefer to avoid storing one of the two files in memory if possible.
EDIT: Here's some pseudocode of exactly what I'm trying to do:
# input: in.bam, in.fastq, chunk_size
i <- 0
for fastq_read in in.fastq:
    bam_read <- extract fastq_read.read_name from in.bam
    n <- i modulo chunk_size
    write fastq_read to out.n.fastq
    write bam_read to out.n.bam
    i <- i + 1

I could swap the above to iterate through the bam file and fetch from the fastq if that's easier.

Does your BAM file already happen to have the same ordering as the fastq files (some aligners produce that by default)? If that's not the case you're pretty much left with (1) sorting files or (2) hashing read names.

Unfortunately they're not in the same order, and additionally some reads exist multiple times in the BAM (supplementary alignments) and a handful of reads do not exist in the BAM at all.

Yikes, your only option is something like what Michael Hall posted then.

Hmm, it's hard to think of a super efficient way of doing this (assuming the files aren't ordered the same - if they are then this whole answer is basically redundant). And also assuming the read ids for both files aren't a perfect intersection.
Off the top of my head you probably want to build a set of read ids for the fastq file and another for the bam. Some python code to get started on that:
import pysam
import itertools

def get_read_id_fastq(ref_path):
    """Extracts the read ids from a fastq file."""
    read_ids = set()
    with open(ref_path, 'r') as ref:
        for line in ref:
            if line.startswith('@'):  # i.e if line is header
                # split the line on spaces, take the first element, remove @
                read_id = line.split(' ')[0].replace('@', '')
                read_ids.add(read_id)

    return read_ids

def get_read_id_bam(ref_path):
    """Extract the read ids from a BAM file."""

    read_ids = set()
    with pysam.AlignmentFile(ref_path, 'r', ignore_truncation=True) as ref:
        for read in ref:
            # query_name is the query template name
            read_ids.add(read.query_name)

    return read_ids

fastq_ids = get_read_id_fastq(fastq_path)
bam_ids = get_read_id_bam(bam_path)

Then take the intersection of these two sets and that's your common read ids.
common_ids = fastq_ids & bam_ids

The next part will be a bit more involved. You will have to iterate through each file one at a time. I would suggest that for the first one you iterate through, create a running dictionary with the read id that was written as the key and the 'chuck number' it was written to as the value. You could create a cycle for your chunk size to manage this effectively
chunk_cycle = itertools.chain(*zip(range(chunk_size)*len(common_ids)))
write_idx = {}

The next part will probably require you to put in some hard-fought times getting it to work. I'll put in some rough pseudo-code to give you an idea.
for line in file:
    read_id = # get line read_id
    if read_id in common_ids:
        chunk_num = chunk_cycle.next()
        write_idx[read_id] = chunk_num
        file_to_write_to = 'out.{}.bam'.format(chunk_num)
        # open this file or write to it if already open

When you go to do the next file when you find a read id is in the common set, you look up it's value in write_idx and this will give you the chunk number to write to.
The reads wont be in the exact same order between the two files, but you could sort later if you needed it (not sure it would matter?).
Hope this helps. Sorry I couldn't give you more, but this should give you a head start hopefully.

Thanks! In fact, I wouldn't need the first part, since every read in the BAM came from the FASTQ. If there are FASTQ records that don't exist in the BAM, I don't mind.

Update: anyone interested can grab my implementation of this solution at https://github.com/scottgigante/nanopore-scripts/blob/master/split_bam_and_fasta.py


text:
What methods exist to calculate RNA expression profile similarity
Some of the work in our lab requires a comparison of a strain across several experimental conditions. We are looking to identify most similar experimental conditions based on the gene transcription response similarity from the cell.
While we could easily invent and create home-grown methods to do it, their implementation and testing are a laborious project in themselves and are outside the scope of our current work.
Are there any methods for RNA expression profile similarity computation that has been already published? If yes, what is your experience using them?

The suggestion to use proportionality is probably the correct one if you are interested in similar patterns between samples. But not if you are interested in absolute differences.
For example: the following two samples are similar in pattern, but not similar on absolute levels:
        Sample 1    Sample 2    Sample 3
Gene A    10          100           80
Gene B     8           80          100
Gene C    12          120          120

Samples 1 and 2 have a perfect proportionality (phi is 0) and also a perfect correlation (as a side note any pair with a perfect proportionality will always have a perfect correlation). However, in terms of logfold changes, samples 2 and 3 are more like each other.
Of course in real life you would never see a comparison like the sample 1 - sample 2 one because normalisation would have removed the scale difference. This was exactly the point brought up by the proportionality paper. But normalisation methods don't normally guarantee that the sum of expression for each sample is identical, and such differences can still occur.
An alternative that might suit more in the second case is either the euclidean distance between the samples, or the euclidean distance on the first two components of a principle component or multi-dimensional scaling. The latter is effectively using mean logFC between samples.
Assuming that x is a matrix containing normalised, log transformed expression values, you could use R and limma to calculate distance in multi-dimensional scaled space as follows:
library(limma)
mds <- plotMDS(x, plot=FALSE)
mds <- data.frame(mds$x, mds$y)
distances <- dist(mds)

Normalisation would only remove the scale difference if it were consistent across a large proportion of genes. In the example you have given, this would indeed be the case if only the indicated genes were assayed, but with this low number of genes it would be more common to normalise based on a housekeeping gene rather than the total gene expression.

Yes, this was my point - hence why there might be reason to use some distance measure (such as euclidean distance) that measured scale difference rather than pattern differences.

There's a new paper that's just appeared on the subject of "proportionality", including a method by which RNA expression might be compared.
This is a new concept for me, and the article is not easy enough for me to read that I can write up a quick summary; the authors don't seem to devote a section in the paper to defining "proportionality". However, here's an interesting chunk from the article:

We graphed the network of relationships between these mRNAs (S5 Fig.), an approach similar to gene co-expression network [12] or weighted gene co-expression analysis [13] but founded on proportionality and therefore valid for relative data. The network revealed one cluster of 96, and many other smaller clusters of mRNAs behaving proportionally across conditions.
  ...
  We are also keen to raise awareness that correlation (and other statistical methods that assume measurements come from real coordinate space) should not be applied to relative abundances. This is highly relevant to gene coexpression networks [12]. Correlation is at the heart of methods like Weighted Gene Co-expression Network Analysis [13] and heatmap visualization [14]. These methods are potentially misleading if applied to relative data.

proportionality is defined as (1+ beta^2 - 2betaabs(r))
where r is the Pearson correlation of log x and log y, and beta^2 is var(log(y))/var(log(x)).
The authors also note: "Proportionality is appropriate, but ϕ does not satisfy the properties of a distance—most obviously, it is not symmetric unless β = 1". They go on to say "Hence, our approach to forming a dissimilarity matrix is simply to work with ϕ(log xi, log xj) where i < j, in effect, the lower triangle of the matrix of ϕ values between all pairs of components."

But why is it called "proportionality" (in non-formulaic terms)? What is the effect of that formula on expression values? Can you fill in the remainder of this sentence: "Proportionality is a measure that represents how the expression of two ..."

No, not really - I didn't quite understand the use of proportionality. Although I do note that phi is not a measure of proportionality, but a measure of goodness-of-fit to proportionality.

Proportionalty implies the ratio of expression values of two genes will be a constant. Instead of fitting $y=mx+c$ in regression, proportionalty fits $y=mx$

If you have a lot of samples with very different environmental conditions, then a Weighted Gene Correlation Network Analysis (WGCNA) might be appropriate.
This type of analysis looks for genes that trace similar (or opposing) expression patterns throughout different conditions (e.g. high-medium-medium-low-absent-high would be highly negatively correlated with low-medium-medium-high-very high-low).
That particular paper introduces the concept of "modules", which are groups of genes that share similar expression patterns. Functions are available for plotting how the expression of canonical module members changes throughout the different conditions, and for identifying which module (or modules) a particular gene is likely to be a member of.
WGCNA works best when there's a lot of different expression changes in the different conditions, which sounds like it would fit well with your project. However, it concentrates more on the genes, rather than the conditions (which seems like it would be less useful for you).

The euclidean distance is probably the simplest, both conceptually and in terms of implementation. It's by no means an elegant solution, and it may not perform well in certain circumstances.
Euclidean distance is easiest to conceptualize as the distance between two points in a two dimensional space.
Y
^
|
|
|
|            * p = (3, 3)
|
|
|
|        * q = (2, 1)
|
----------------------------> X

In this example, the distance between the two points is
d(p, q) = sqrt( (p_x-q_x)^2 + (p_y-q_y)^2 )
        = sqrt( (2-3)^2 + (1-3)^2 )
        = sqrt(5)
        ≈ 2.24

For a gene expression profile with two genes, this is exactly how the euclidean distance would be calculated, using expression values from one gene as the X axis and expression values from the other gene as the Y axis. Realistically, though, gene expression profiles typically contain thousands or tens of thousands of genes, so instead we use the generalization of the distance calculation for N dimensions.
d(p, q) = sqrt( (p_1-q_1)^2 + (p_2-q_2)^2 + ... + (p_N-q_N)^2 )

Packages for R and Python make these types of calculations trivial once you have the data loaded into the correct data structure. See Ian's answer for some example R code.


text:
What sequencing data metrics should I record?
I'm new to bioinformatics and Im starting a new microbial sequencing project and I'd like to record all of the qc data correctly. This is research and I'm not sure what I'll need later. Is there any consensus on useful qc metrics to record in a database for each sequencing run?
No rnaseq, running FastQC, assembling with a defined reference for a single species.

The metrics from FastQC are more of a snapshot. It depends on your sequencing platform used. Was it illumina, iontorrent, pacbio?

The sequencing platform is illumina. Primarily miseq. I just don't want to regret not recording information that might be useful later.

Well all the important info for illumina will be stored in basespace you would be looking at cluster density and the Q-score.

Is there any reason not to keep all of it?

No, just wanted to distill the essentials into a local database if possible.

What it might not be essential today might be essential in a couple of years. You can always delete it later, but you can't recover data after years of the experiment

If you are using a MiSeq is your data stored on basespace or locally?  If it is on basespace then the QC metrics should be store there for you.  If locally then you can go into the run data and pull the metrics files from the Illumina reports.  Really you are looking for information about clustering etc.  The sequencing quality and read metrics can be recalculated from the fastq/Bam files at any time.

@AaronBerlin we're storing locally. Thanks for the suggestions.

This also depends what you want to do with the data. The FASTQs already have all the quality scores stored, so it's easy to re-create all the QC outputs if necessary. MultiQC is good at making these reports: http://multiqc.info/
If you're planning to assemble them, the metrics matter a bit more. You'd want to record things like coverage depth, number of contigs, length of contigs etc. People usually report N50 for example

I would say that the following metrics are useful:

Sample source (Extraction place, person, date)
Details about pre-HTS preparation (e.g. amplicon PCR, rRNA removal)
Reads per sample (+ per barcode / lane, if available)
Sequencing kit (e.g. TruSeq v3), useful for adapter trimming
Sequence length (e.g. 125bp)
Paired end / single end
If paired end, target fragment length in sample prep (e.g. 650bp)
Reference sequence metadata

Scientific name
Location of reference genome (preferably a stable web address)
Version of reference genome
Read mapping rate (per sample / barcode / lane), as a percentage of total reads
Mapping program + version

If available, potential contaminating organisms

Name / description
Location of contaminant sequence (preferably a stable web address)
Estimated mapping rate (bearing in mind that there may be some overlap with the target)
Mapping program + version

More generally, ENA / EBI has an excellent reference in sample checklists, which indicate the minimum acceptable metadata for attaching to sequence archive data:
https://www.ebi.ac.uk/ena/submit/checklists


text:
Alignment based vs reference-free (transcriptome analysis)?
I want to focus on transcriptome analysis. We know it's possible to analyze RNA-Seq experiment based on alignment or k-mers.
Possible alignment workflow:

Align sequence reads with TopHat2
Quantify the gene expression with Cufflinks

Possible reference-free workflow:

Quantify sequence reads with Kallisto reference-free index

Both strategy generate gene expression table.
Q: What are pros and cons for each of the approach? Can you give guideline?

Reference-free quantification doesn’t really exist (and Kallisto isn’t it).

Here is a project that tries to build tools for high throughput sequencing data without a reference: https://colibread.inria.fr/project/

First of all, I would emphasize that "alignment-free" quantification tools like Salmon and Kallisto are not reference-free. The basic difference between them and more traditional aligners is that they do not report a specific position (either in a genome or transcriptome) to which a read maps. However, their overall purpose is still to quantify the expression levels (or differences) of a known set of transcripts; hence, they require a reference (which could be arbitrarily defined).
The most important criterion for deciding which approach to use (and this is true of almost everything in genomics) is exactly what question you would like to answer. If you are primarily interested in quantifying and comparing expression of mature mRNA from known transcripts, then a transcriptome-based alignment may be fastest and best. However, you may miss potentially interesting features outside of those known transcripts, such as new isoforms, non-coding RNAs, or information about pre-mRNA levels, which can often be gleaned from intronic reads (see the EISA method).
This paper also has some good considerations about which tools may work best depending on the question you want to answer.
Finally, another fast and flexible aligner (which can be used with or without a reference transcriptome) is STAR.

I wouldn't say Kallisto (or Salmon) are reference-free. They use a transcriptome as reference anda concept called pseudo-alignment which greatly speed up the process of assigning your reads to a transcript.
That said, both approaches of (i) mapping against a reference genome (what you called alignment workflow ) and (ii) mapping against a reference transcriptome will serve different purposes
Transcriptome mapping using pseudoalignent is becoming the method of choice for gene/transcript quantification and differential expression analysis. The drawback is that you only focus on known transcripts
Two typical workflows are:

Kallisto followed by sleuth
Salmon, followed by tximport and DESeq2/EdgeR

Genome mapping is useful for, per example discovery of new isoforms. You shouldn't use TopHat anymore as it has been discontinued by the author.
A typical workflow would be:

Hisat2 (alignment)
StringTie (transcript assembly and abundance estimation)
Ballgown (differential expression)

It's just worth noting (in addition to this very useful answer) that Salmon can also be used with sleuth by means of the wasabi package.


text:
How to build a BAM header file with htslib in C++?
I'd like to use C++ to generate a new BAM file programmatically. This is an example how to use htslib to generate a new BCF file on the fly.

https://github.com/samtools/htslib/blob/develop/test/test-bcf-translate.c

The two key functions are:

bcf_hdr_append
bcf_hdr_add_sample

Now when I navigate to the header file for SAM/BAM functions:

https://github.com/samtools/htslib/blob/develop/htslib/sam.h

The only "append" function is bam_aux_append. But that's for optional fields.
My questions:

Given a bam_hdr_t header (potentially empty), how to add a new chromosome reference? (For example, add a new chrAB to the BAM file)
Given a non-empty bam_hdr_t * header, how to remove a chromosome? (For example, remove chr4 from the BAM header)

My sample code:
// Empty BAM header, nothing inside
bam_hdr_t * hdr = bam_hdr_init();

// I'd like to add a new synthetic `chrAB` into the header, how?

You're just looking in the wrong header file. The BCF stuff will be in vcf.h (you'll see the two functions you mentioned there).
With BAM files one doesn't typically need to add/remove headers, which is why there are not convenience functions for them. To add a header:

Increment hdr->ntargets
realloc(hdr->target_len, ...) and hdr->target_name.
You likely want to append the name/length to the end of hdr->target_name and hdr->target_len, respectively.

I don't think you need to update the dictionary or plain text for most uses (or at least I've never needed to).
Removing a chromosome would follow the reverse instructions, (you'd need to shift everything in the name and length arrays if you remove entries in the middle.
WARNING, doing this is prone to error!
I would strongly cautious you to be very careful when doing this. Please remember that the chromosome that entries are mapped to isn't actually stored in alignment entries themselves. Rather, they simply have an index into hdr->target_name, which saves space. Consequently, if you start adding/removing headers to/from a BAM file with alignments then you're going to have to modify all of the alignments to keep the correct alignment:chromosome associations. Otherwise you'll end up with completely useless results (this is the same caveat as with samtools reheader).


text:
Is this a low complexity region in our human genome?
I have a screenshot where many of my reads were aligned to a region that I suspect a low complexity region. Although you can't see, all those reads are clipped in the cigar strings. Sample cigar strings: 64S32M29S, 74S32M18S etc... Consequently, the actual sequence mapped to the genome is less than the read length.
I have a feeling that my alignments are bogus because of the complexity and clipping, but I'm not sure if the complexity is indeed low. All I see is a few "T" between a bunch of "A". How to define a low complexity region?
Q: Is this a low complexity region in the human hg38 genome? Would that be fair if I report "the alignments are likely be errors due to the low complexity and alignment clipping"?

Yes, that's a short low complexity region wedged between a SINE and an snRNA. More importantly, your alignments have MAPQ of 0 (that's why they're filled with white in IGV), which will happen if they map equally well to multiple locations. Without looking at the sequence one can use that alone to determine that these are not trustworthy mappings.

Oh, I didn’t know about white meaning MAPQ=0. Very useful!

Yes, this is a low complexity region.
Regions are considered low complexity (or having a simple sequence) when they contain an abundance of a single base, or an abundance of short tandem repeats. The simplest is a tandem repeat of a single base (e.g. AAAAAAAAAAAAAAAA), also called a homopolymer. In your case there is a 3-unit tandem repeat of TAAAAA, and surrounding As (with a few intermingled Ts). Depending on the genetic history of this region, it might be that those other As and Ts are derived from the same TAAAAA unit, with some slipping that has occurred during replication to create a cryptic (i.e. non-perfect) repeat structure.


text:
Read sorted and indexed BAM files faster in C++?
I have some sorted and indexed alignment large BAM files. I'd like to read all the reads, and apply custom operation on each of them. The orders aren't important. I'm using htslib in C++ to do the reading. I have a single powerful machine.
The problem is ... reading a large BAM file is very slow even in C++.
Q: Is there any trick to speed up the reading? Reading BAM files multithreaded (I have many cores)?

If you're reading the entire file, indexing won't help. That's for quickly seeking to specific regions in a BAM file.

@gringer Can I speed up by splitting the file into partitions, then each of those by multithreading? I have many cores to do that.

Reading data is an I/O problem, not a CPU one. Do you have your data on an SSD?

@AlexReynolds Yes. But is there anything else we can do in programming?

Using mmap routines to read data, maybe? It might give you a slight benefit over reading with fopen/fseek/fread etc. But there's no way to guarantee a benefit, short of benchmarking whatever it is you are trying to do. And you have to be careful about purging caches between tests. It's easy for people to make spurious claims about their tools being fast if they aren't purging the read cache.

You might investigate terms like "striping" or "parallel I/O" to gain benefit from multiple read operations via hardware. Having multiple processors or cores can be useful for algorithms, but I don't think it will benefit reading sequential data from storage.

The optimal solution is entirely determined by the "custom operation" you want to apply. What do you want to do exactly?

What do you mean by "very slow"? How does it compare with time cat <file.bam> > /dev/null? How does it compare to the code when nothing is done with the records?

In this answer I want to show you some benchmarks that compares three different serial ways of reading data in C++, the third one being the fastest.
In the first method I use an std::istreambuf_iterator, in the second I read the file line by line into a std::vector and in the third one I use C flavoured operations.
/*
compile with g++ -std=c++11 -O2 reading_files.cpp -o reading_files
test the program with the following
dd if=/dev/random of=./dummy5GiB bs=1024 count=$[1024*1024*5]
 for i in 0 1 2; do time ./reading_files dummy5GiB $i; done
*/

#include <fstream>
#include <iostream>
#include <vector>
#include <cstring> //for strcmp
#include <string>
#include <algorithm> //for min

using namespace std;

int main(int argc, char * argv[])
{
    if(argc != 3) {
        cerr << "usage: reading_files <path> <number 0 or 1 or 2>\n";
        return 1;
    }

    ifstream file(argv[1], ios::in|ios::binary);
    if(!file.is_open()) {
        cerr << argv[1] << " file not opened!\n";
        return 1;
    }
    cout << "---------------------------\n";
    if(!strcmp(argv[2],"0")) {
        cout << "method 0\n";
        vector<char> content((std::istreambuf_iterator<char>(file)), std::istreambuf_iterator<char>());
        cout << content.size() << " bytes read\n";
    } else if(!strcmp(argv[2],"1")) {
        cout << "method 1\n";
        string line;
        vector<string> content;
        while(getline(file,line)) {
            content.push_back(line);
        }
        cout << content.size() << " lines read\n";
    } else if(!strcmp(argv[2],"2")) {
        cout << "method 2\n";
        file.seekg(0, std::ios::end);
        unsigned long long file_size = file.tellg();
        cout << "file size is " << file_size << "\n";
        //maximum 5 GiB of memory, so to fit my RAM
        unsigned long long max_block_size = 1024ULL * 1024ULL * 1024ULL * 5;
        file.clear();
        file.seekg(0, std::ios::beg);

        unsigned long long block_size = min(max_block_size,file_size);
        char * buffer = (char*)malloc(block_size*sizeof(char));
        unsigned long long processed = 0;

        while(!file.eof()) {
            file.read(buffer, block_size);
            cout << "read a block of " << min(file_size - processed, block_size) << " bytes\n";
            //do work here with the block in memory (if the block size is > 0)
            processed += block_size;
        }

        free(buffer);
    }
    file.close();
    return 0;
}

In the third method I read the file in blocks of 5GiB, so to fit the free RAM of my machine.
Here is the benchmark of the program on a 5 GiB randomly generated file.
bash-3.2$ for i in 0 1 2; do time ./reading_files dummy5GiB $i; done
---------------------------
method 0
5368709120 bytes read

real    0m35.194s
user    0m19.665s
sys 0m12.999s
---------------------------
method 1
20973135 lines read

real    0m40.478s
user    0m34.046s
sys 0m5.870s
---------------------------
method 2
file size is 5368709120
read a block of 5368709120 bytes
read a block of 0 bytes

real    0m4.413s
user    0m1.757s
sys 0m2.617s

Here is the benchmark on a real 32GiB .bam file using only the third method, here the splitting of the file read plays a fundamental role. I can't tell you how big because my hard disk memory gets saturated before the end of the execution and so I have to kill the program.
time ./reading_files ./HG00252.mapped.ILLUMINA.bwa.GBR.low_coverage.20130415.bam 2
---------------------------
method 2
file size is 34316054058
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 2103799338 bytes

real    0m50.656s
user    0m9.135s
sys 0m35.728s

The key point is to minimise the file operations reading large block of data at once (as long as they fit the memory) and operate on them.
Using mmap could lead to an improvement, but I don't think it would be as large as using the third method instead of the first two.

Why is your “method 2” using gratuitous manual (C-style!) memory management instead of a vector? — I’m also not sure how much this actually helps OP since the question is about parsing BAM, and as far as I know neither Staden IO nor htslib supports reading from a custom
buffer (but should perform its own buffering anyway).

If you need to iterate through the entire BAM file, decompression is usually the main bottleneck.  Three things you can do:

Install libdeflate (https://github.com/ebiggers/libdeflate ) and then the development version of htslib, and make sure the latter is compiled to use libdeflate.  This roughly doubles decompression efficiency in my experience.
Call htslib's bgzf_mt() function right after opening the file for reading, to request multiple decompressor threads.  I generally go with 2-6, depending on what I'm actually doing with the reads.
Have a separate read/decompress-ahead thread perform the bam_read1() calls.  Ideally, when your main thread is processing read #X, background threads are already working on decompressing reads e.g. (X+100) through (X+200).

Also, if your workload is "embarrassingly parallel", you may also benefit from launching multiple processes simultaneously, configuring each to handle a similarly-sized chunk of the BAM.

I stumbled upon this when looking for ways to improve my BAM processing and I was wondering if you have a code example that you could share?

In both deepTools and MethylDackel we implement multithreading by having individual threads process separate genomic regions. At some point you hit an IO limit, but unless your processing is trivial then IO won't typically be limiting with single-threaded applications. Since you already have sorted/index BAM files, you can take advantage of the same strategy.


text:
How to filter out cross alignments from a BED file?
I have a BAM file:
@SQ SN:chr1 LN:248956422
@SQ SN:chrx LN:248956423
ST-E00110:348:HGVKKALXX:1:1201:5822:48670   323 chr1    9999    0   67H66M16H   chrx    1000    0   GATAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC  JJJJJJJJJJJJJJJJAJJJJJJJJJJJJFJJJJJJFJFJJJJJJFJJJJJJJJJJJA77FJFJJJ  NM:i:0  MD:Z:66 AS:i:66 XS:i:65 SA:Z:chr5,18606834,-,73S76M,34,0;   RG:Z:g1

There is a read aligned to chr1, and it's mate aligned to chrx.
I have a BED file:
chr1    0   100000  TestOnly

I would like to filter out everything that falls outside my BED region, that includes cross-alignments. In my example, although my read aligned to chr1 but it's mate is not. I don't want this read.
When I do:

samtools view -L test.bed test.bam

the command gives me the read because it doesn't check cross-alignments.
My solution:

samtools view -L test.bed test.bam | grep -v chrx

but this is very slow and clumsy. In my production pipeline I would have to do something like:

samtools view -L test.bed test.bam | grep -v chrx | grep -v ... | grep -v ... | grep -v ... | grep -v ...

Q: Is there a better solution?

According to the SAM specification, the 3rd field of a SAM line (RNAME) is:

RNAME: Reference sequence NAME of the alignment. If @SQ header lines
  are present, RNAME (if not ‘*’) must be present in one of the SQ-SN
  tag. An unmapped segment without coordinate has a ‘*’ at this field.
  However, an unmapped segment may also have an ordinary coordinate such
  that it can be placed at a desired position after sorting. If RNAME is
  ‘*’, no assumptions can be made about POS and CIGAR.

And the 7th field is (emphasis mine, missing "to" theirs):

RNEXT: Reference sequence name of the primary alignment of the NEXT
  read in the template. For the last read, the next read is the first
  read in the template. If @SQ header lines are present, RNEXT (if not
  ‘*’ or ‘=’) must be present in one of the SQ-SN tag. This field is set
  as ‘*’ when the information is unavailable, and set as ‘=’ if RNEXT is
  identical RNAME. If not ‘=’ and the next read in the template has one
  primary mapping (see also bit 0x100 in FLAG), this field is identical
  to RNAME at the primary line of the next read. If RNEXT is ‘*’, no
  assumptions can be made on PNEXT and bit 0x20

So, you want to remove those lines whose 7th field isn't = and, just in case, those lines whose 7th field isn't = and isn't the same as the 3rd field. You can therefore use something like this:
samtools view -L test.bed test.bam | awk '$7=="=" || $3==$7

And, to save as a bam file again:
samtools view -L test.bed test.bam | awk '$7=="=" && $3==$7 |
    samtolls view -b > fixed.bam

On a separate note, there's very rarely a need to chain multiple grep commands like that. You can just use \| (or | with the -E or -P options) to separate them. Something like:
samtools view -L test.bed test.bam | grep -v 'chrx\|chr2\|chr10\|chrN'

Or
samtools view -L test.bed test.bam | grep -Ev 'chrx|chr2|chr10|chrN'

If you do it like this, the fixed.bam file is missing the header, which in my experience generates a lot of problems. I recommend always adding the header back in; either by specifying -h when reading the original BAM, or by adding it separately: (samtools view -H infile.bam; samtools view …) > samtools view -b > outfile.bam.


text:
How to filter allele frequency VEP with bcftools?
I have a VEP annotated VCF file, and I wish to use bcftools to filter it. bcftools work very well on standard VCF fields, but I'm not able to make it work on VEP fields.
For example, the following command would report "not defined in the vcf header" (it's not a standard VCF field but part of the VEP fields in the VCF file). split-vep is a VEP plugin offered by bcftools.
bcftools +split-vep MyFile.vcf -i 'gnomAD_exomes_AF<0.01'

Thanks for the edit! Now, please also give us a file we can test on. Try to make it as minimal as possible while reproducing the error. Note that if the field isn't actually included in the headers, then the VCF is not valid, and bcftools won't know how to deal with it. You might, however, be able to script something, but for that we will need to see at least some part of the file.


text:
How to interpret duplication from MultiQC/FastQC?

FastQC of my sample files, aggregated into a single plot by MultiQC. Blue represents unique reads. Black represents duplicate reads. The x-axis is the number of reads.
I see I have more duplicate reads than unique reads, how do I interpret this plot? This is my attempt:

FastQC took the first 100,000 sequences in each sample (details).
  Generally, higher duplicated read percentage is welcome as it
  represents better sequencing depth and lower sequencing error
  contamination. This dsitribution in this plot looks reasonable and
  nothing unusual.

Can you provide us with some more info, including the typical number of reads in each library, the library type, source material, et cetera? It's hard to tell how to interpret your data without that info. If it's single-cell versus a sample from bulk tissue that makes a big difference in complexity. Also, something like RNA seq will have more duplicates that genome sequencing library types.

I would agree with your statement.  Typically the ratio of unique to duplicated reads is a function of input library complexity (number of input molecules) verse sequencing depth (the amount of molecules sequenced).
In the example you gave you have about a 2:1 ratio of duplicates to unique reads which suggests that these libraries have been sequenced to their fullest and additional sequencing would yield few new unique reads.

The higher number of duplicates could be in a high-complexity library sequenced very deep or in a low-complexity library sequenced with many fewer reads. Without more info from OP it is hard to interpret.

the x-axis provides the total number of reads.  either way it is the ratio that is important to tell you if the sequencing has saturated the library.


text:
How should I address batch effects in my experiment?
Let's say I have an RNA-Seq experiment, where I'm interested in the significantly differentiated genes between pre-treatment and post-treatment conditions. "rep" == biological replicate.
Sample PreA (3 reps **all** in 2017) vs Sample PostA (3 reps **all** in 2018)
Sample PreB (3 reps **all** in 2017) vs Sample PostB (3 reps **all** in 2018)
Sample PreC (3 reps **all** in 2017) vs Sample PostC (3 reps **all** in 2018)

I also have a sample:
Sample PreD (3 reps all in 2017) and Sample PreD (3 reps again in 2018)

It is safe to assume the samples in sequenced in 2017 were done differently to 2018.
Q1: Is there problem with batch effects in my study? When I do a test, there should't be any way to split variation from batch effects and biological variation?
Q2: If Q1 is correct. How should I reduce batch effects?

Can I add a batch factor to my model in DESeq2/edgeR?
Should I use my PreD control samples to derive normalisation factor as in RUV? (assume PreD as a negative control)
Should I use sva package?

EDIT
Devon suggested for adding a batch variable, but how? That's how I'm doing (ignoring batch effects):
dds <- DESeqDataSetFromTximport(txi, data, ~Sample)

data holds a data frame with a Sample column. Sample would be "PreA" and "PostA" in the first test. Very simple.
I couldn't add Date to it, like:
dds <- DESeqDataSetFromTximport(txi, data, ~Date + Sample)

because the code crashed complaining of the 2017/2018 problem discussed above. I'm not sure how I should add PreD to the design.

Yes, though thankfully Your 2018 PreD samples will help you resolve this.
Simply add the batch effect to the design (~Batch + Treatment) and DESeq2 (or edgeR or Limma) will handle this for you.

You do not need SVA or RUV, thankfully, since you quite cleverly sequenced one group in both batches.
To clarify, your coldata will be something like:
Group    Time    Batch
A        Pre     A
A        Pre     A
A        Pre     A
A        Post    B
A        Post    B
A        Post    B
...
D        Pre     A
D        Pre     A
D        Pre     A
D        Pre     B
D        Pre     B
D        Pre     B

And the design is actually ~Batch + Group + Time. You can probably add an interaction term in if you then modify the model matrix manually to ensure it's full rank.

Thanks. PreA vs PostA, involve no "PreD" at all. How should I address it? "PreD" is a sample name, not like a factor like as.factor(c("2017", "2018")).

You're really interested in the A:Post interaction. If you phrase it and think of it that way the model matrix is easier to construct (see my update to get you started).


text:
how to find the bound form of an enzyme structure?
For my undergraduate research I'm looking for a database that gives the bound form of a particular protein structure. Is there any database that provide us with such data? So far I've found following proteins related to my work

Ascorbate peroxidase (APX)
Beta-glucosidase
Calcineurin b-like protein-interacting protein kinase.

I've searched in PDB and UNIPROT.

Hi Lock23, thanks for your question and welcome to Bioinformatics Stack Exchange. After reading your question, I did a quick check on the Biology and Chemistry Stack Exchanges, and was surprised at how little "database" stuff there was.

For experimentally determined structures PDB is the database and it's unlikely you'll find structures that are publicly available but not in the PDB. There is also protein structure prediction software (servers).

Tried looking for an explicit database? i.e.

ComSin: database of protein structures in bound (complex) and unbound (single) states in relation to their intrinsic disorder: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2808974/ or
LigASite: a database of biologically relevant binding sites in proteins with known apo-structures https://www.ncbi.nlm.nih.gov/pubmed/17933762

Else the PDB itself has a very powerful advanced search that should allow you to do this: http://www.rcsb.org/pdb/search/advSearch.do

Are you aware of BRENDA? I was just introduced to it today for a completely separate reason (looking at carbohydrate enzyme families in the Nippostrongylus brasiliensis proteome), and it seems to be a fairly comprehensive database. There is at least a literature link there for the ascorbate-complexed crystalisation of Ascorbate peroxidase on that site. Following through the paper, I see that the complex has been uploaded to PDB.


text:
predict the foldability of single-stranded DNA molecules
I have a list of regions of the human genome and I want to predict if single-stranded molecules in a buffer would tend to fold and create pin structures by sequence self-complementarity. What's the most precise software/parameter set I can use to such effect?
TRIED SO FAR:
Vienna RNALFold (unsure if applying an RNA software to DNA is adequate):
Compute locally stable RNA secondary structure with a maximal base pair span.
For a sequence of length n and a base pair span of L the algorithm uses only
O(n+L*L) memory and O(n*L*L) CPU time. Thus it is practical to "scan" very
large genomes for short RNA structures.
Output consists of a list of secondary structure components of size <= L, one
entry per line. Each output line contains the predicted local structure its
energy in kcal/mol and the starting position of the local structure.

Example:
samtools faidx $ref $region | seqtk seq -A - | $HOME/viennarna/bin/RNALfold --noconv
>chr1:100000-100150
.((((.(......).)))). ( -2.70)  127
.(((.((........)).))). ( -5.70)  123
.(((((((.((........)).)))..)))). ( -6.30)  119
.(((.(.(((.((.....))))).)..))). ( -7.00)  117
.((((.((........)).)))). ( -8.30)  109
.((((..((.((..(..(((.....)))..)..)))))))). (-10.10)  103
.((((.((..((..(((...((((.((........)).)))).....)))))..)))))). (-10.50)   90
.((.((.((((....................)))).)).)). ( -4.10)   77
.((((.((.((((....................)))).)).)))). ( -5.70)   75
.((((..((((.((....)).)))).)))). ( -2.40)   72
.((.((..((((..((((.((....)).)))).))))..)))). ( -5.50)   65
.(.((..(..................((.((..((((..((((.((....)).)))).))))..))))..((((.((........)).)))).)..)).). (-16.20)   40
.((((......(.((..(..................((.((..((((..((((.((....)).)))).))))..))))..((((.((........)).)))).)..)).).....)))). (-18.10)   30
.(((....(((((((((...(.........)..)))).)))))((.((..((((..((((.((....)).)))).))))..))))..((((.((........)).))))......))). (-18.50)   23
.((....(((....(((((((((...(.........)..)))).)))))((.((..((((..((((.((....)).)))).))))..))))..((((.((........)).))))......)))....)). (-18.70)   17
.(((((((((......((...(((((..((((...((...))....))))..)))))...))......))))))))). ( -9.60)    7
.((.(((((((((......((...(((((..((((...((...))....))))..)))))...))......))))))))).)). (-10.90)    4
.(((..(((((........)))))..))). ( -3.60)    1
ACTTAAGTTGTAGAAGGGAATAACGCAAGAGTGAATTTAGGGCGGGGCAAAAGGATAAATTTTACGGTACAAAGTTTCTACGGGTTTATATATGTATAACTAAGTCCAAGCGCGGGGGATATGGCCAGTGCACAACGGCGGGCATCATAAT
 (-21.70)

I think this is a nice question but a bit open-ended. Are you searching for other tools because this didn't convince you ? Or what requirements are you looking in a other tools? How do you define "best"? Less FDR or more precision ?

I edited my question to "most precise" and added that the reason I am unsure about RNALFold is because it's an RNA software applied here to DNA molecules.

Great edit, I didn't realize you were using DNA as input, I thought you were using RNA from those regions

You could also take a look at Infernal.
I think your best bet is to use RNA-focused programs; I'm not aware of a DNA-specific one, which makes some sense because DNA is typically double-stranded, so the secondary structure community is primarily focused on RNA. I would predict that secondary structure potential of single-stranded DNA would be similar to that of RNA.

Thanks for recommending. I learned about Infernal about 5-8 years ago, and after reading the documentation as it is now, it seems to me that there isn't an easy way of asking for the folding of DNA sequences. The commands seem to be based on building alignments. Am I wrong?


text:
remaining human genome variation that hasn't been sequenced?
Given all the genome variation information we have accumulated from resequencing human genomes up until now (early 2018), how much variation is left that hasn't been sequenced?
EDIT: given one of the answers, and the fact that 99% of the SNPs at 1% or more are found, is there a quantitative measure of remaining variation?

I don't think this is answerable. Do you mean single nucleotide variants? Copy number changes? In healthy individuals only? Only living individuals? All possible variants? Pretty much anything is possible. . .

Please clarify "variation". If you limit this to single nucleotide polymorphisms only and to variants above a certain frequency threshold, as the answer you have received has done, then it might indeed be answerable.

In theory, almost any base in the human genome may mutate, so you have billions of variants to go. Ok, this is not so useful.
A related and potentially useful question is: given a human, what is the fraction of his/her variants seen in the the previously sequenced samples? For a Wright-Fisher population, there is an analytical answer: x% of variants on a haplotype has frequency below x%. Replace "x" with "1": 1% of variants on a haplotype have frequency below 1%. The 1000 genomes project concluded that 99% of SNPs at frequency 1% or higher have been called by the project. If we all came from a Wright-Fisher population, 99% of of SNPs have been called in 1000g. Given that the project also calls SNPs at lower frequency, nearly all SNPs in a newly sequenced sample are present in 1000g.
The above is a theoretical analysis. In practice, there are two major complications. First, we are not a Wright-Fisher population. Due to recent population expansion, each haplotype harbors more new mutations than the Wright-Fisher model would predict. Given a known population history, it is actually possible to numerically compute $f(x)$, the fraction of variants with frequency below x (for Wright-Fisher, $f(x)=x$). I don't have this result, though. Second, 1000g and many other genome projects didn't call SNPs in repetitive regions and missed a significant portion of high-frequency indels and structural variations. With Illumina, you can only call ~30% of long deletions callable with PacBio data, but we only have a dozen of PacBio human genomes in public. There are still a lot we don't know.


text:
creating population AF from set of individuals in 1000 genomes population?
What's the recommended way to create population-level Allele Frequency variant entries from a set of individuals (vcfs?) in the 1000 genomes (or similar) project?
For example, given the IBS individuals in the 1000 Genomes project:
http://www.internationalgenome.org/data-portal/population/IBS
I would like to create an IBS .vcf.gz that contains an IBS population-level AF entry for each of the variants present.
Any ideas on how to achieve this?
EDIT: this is related to my earlier question on how to create the most complete file for population allele frequencies, here:
what is the most complete vcf file of population allele frequencies that can be built/downloaded from public datasets?
I have seen that there is an Allele Frequency Calculator tool in EnsEMBL, which is at least available with GRCh37:
http://www.internationalgenome.org/allele-frequency-calculator-documentation
http://grch37.ensembl.org/Homo_sapiens/Tools/AlleleFrequency?db=core
There are also instructions on how to do this genome-wide using vcftools:
vcf-subset -c ~/1000g/IBS.samples.list ~/1000g/ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz | fill-an-ac | bgzip -c > ~/1000g/IBS.chr1.20130502.vcf.gz

This seems to be generating correct files, and the AF can be calculated with a tool that combines AC and AN values, such as it is described here:
https://www.biostars.org/p/123468/
Any other ideas / warnings / pitfalls of the approach?

Are you looking to get this genome-wide? Or are you only interested in particular regions/variants?

Hi @Emily_Ensembl looking to do this genome-wide for as many populations as possible. I see there are 25-30 available in EnsEMBL from the 1000 genomes project.

bcftools has a plugin called fill-tags which should do the job for you easily, without having to mess around with lots of complicated bash scripting.
You could do something along the lines of this, presuming you have a .vcf file with all the IBS individuals in it:
bcftools +fill-tags in.vcf.gz -Oz -o out.vcf.gz

This will give each position in the .vcf file new INFO tags, one of which will be allele-frequency.
You can then pull off the frequency for a particular variant (CHR:POS) by:
bcftools query -t CHR:POS -f'%CHROM\t%POS\t%AF\n' out.vcf.gz

I will check it out. Thanks


text:
Selecting 65000 SNPs where AF is close to 0.5 in all or most populations
I am evaluating the tool somalier (https://github.com/brentp/somalier) and I need to create a list of about 65,000 SNPs where the allele frequency (AF) is as close to 0.5 as possible across the most representative set of populations possible with today's data.
My alignments are in GRCh38, so I need a vcf input file for somalier also in GRCh38, not the one that is available in the somalier website, which is in GRCh37.
I am starting as input with the vcfs from the 1000genomes project in GRCh38 positions, and I would like to come up with a filtering strategy that:
(1) Sorts all SNPs by the AF being close to 0.5 in as many populations as possible.
(2) Takes the first 65,000 SNPs from that list and produces a vcf file from them.
My simplistic strategy so far has been to take the global AF value if it's between 0.3-0.7 and print it if it's a SNP with this perl one-liner:
gunzip -c my.vcf.gz | grep -v '#' | perl -lne 'if ($_ =~ /\;AF=0\.[34567]/) { if ($_ =~ /\w+\t\d+\t\S+\t\w\t\w\t/) {print $_} }'

Any recommendations?

If possible, I strongly urge you to use hg19 instead of GRCh38 because there are much better frequency data available. Basically, GnomAD uses hg19, so any gnomAD frequencies for hg38 will come from lifting over the positions and that doesn't work too well in my experience.

I haven't been able to reliably liftover the sites.vcf.gz file that can be downloaded from the github.com/brentp/somalier site. I will create another Q for the liftover issue.

What do you mean? Somalier's github page says it's for hg19, not hg38: "where sites is a VCF of variant sites (provided by somalier for hg19)."

My alignments are in GRCh38, so I need a vcf input file for somalier also in GRCh38, not the one that is in GRCh37 (gh19).

Just use bcftools view for filtering:
$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz > output.vcf

To truncate this list to 65,000 SNPs count the header lines, sum them to the number of SNPs you like to have and use head -n.
$ bcftools view -h input.vcf.gz|wc -l
255

$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz | head -n 65255 > output.vcf

Using bcftools is the right way, but you could also greatly simplify the command you're using. You can do the same thing with a single zgrep call:
zgrep -P '^[^#]+\t\d+\t\S+\t\w\t\w\t.*;AF=0\.[3-7]' my.vcf.gz

And to limit to the first 65000 results:
zgrep -m 65000 -P '^[^#]+\t\d+\t\S+\t\w\t\w\t.*;AF=0\.[3-7]' my.vcf.gz

Finally, to keep the header of the original vcf unchanged, you could do:
zgrep -m 65000 -P '^(#|[^#]+\t\d+\t\S+\t\w\t\w\t.*;AF=0\.[3-7])' my.vcf.gz

I also suggest you use hg19 instead of hg38 data since the best available source for frequencies, the gnomAD database, is based on hg19 and the tool you want to test also provides a site file for hg19. But that may not be possible for you.
If your VCF file has multiple AF entries for different sub-populations and you want to use those, you'll probably need to script something. We can't really help with that unless you show us some example input lines and start writing the basic script yourself though.

I have a different opinion about that. hg19 is outdated sind 2013. Whenever you start something new, you should move to the current reference genome, which is by now hg38. There are lifted over versions of 1000 Genomes, ExAC and gnomAD which provide at least the frequency data by ensembl. As long as the frequency data are given in the INFO column you can easily filter by them using bcftools. So there is no need to script something.

@finswimmer not sure what you are referring to, but if it's about hg38, I am not saying hg19 is better (it isn't) only that the gnomAD data is based on hg19 and liftover is imperfect. As for the rest, yes of course bcftools is the way to go!


text:
Accessing .bam/.cram files from AWS S3?
What's the best/easiest way for accessing .bam/.cram files from S3?
I have .bam/.bai .cram/.crai and .bed + .gff files that I want to make available easily to others so they can browse in IGV and other solutions from their computers. There seem to be different solutions after cursory Googling, what would you recommend?
Here is what I've found so far:
IGV:

https://twitter.com/igvteam/status/1187766079457644544
https://umccr.org/blog/igv-amazon-backend-setup/

Samtools: samtools directly supports AWS

http://www.htslib.org/doc/htslib-s3-plugin.html

The full URL format is:

    s3[+SCHEME]://[ID[:SECRET[:TOKEN]]@]BUCKET/PATH

samtools directly supports AWS these days: http://www.htslib.org/doc/htslib-s3-plugin.html

You could set up a public read-only bucket in an S3 account and put your files there. URLs would look like:
https://${username}.s3.amazonaws.com/public/${filename}

For samtools, the request for a BAM file's index should be built off that URL. For IGV, I think you can specify separate, custom URLs for the data and the index resources.
Public buckets are available to all, without need for credentials. So that might be an issue to think about, with respect to billing for web requests and data transfer.
I'm still not 100% certain that samtools/tabix support https requests reliably (I have a prior question on here about that), but that might be one approach.

Doo you have your heart set on AWS, or are you open to using Google Cloud Platform (GCP)?
If you're open to using GCP, there are some advantages. For example GATK can read and write directly to and from Google storage buckets, and IGV can use your login credentials to access data in Google storage as described here:
https://googlegenomics.readthedocs.io/en/latest/use_cases/browse_genomic_data/igv.html
Note that if you make the data public, by default you're going to get billed if people download it. Not sure how that's handled on AWS, but on GCP you can set the bucket (or just the file) to be "requester-pays", ie the person who downloads is the person who pays for it:
https://cloud.google.com/storage/docs/requester-pays
In general there are several ways to interact with Google buckets.

web-based file browser (example is a public bucket I manage): https://console.cloud.google.com/storage/browser/genomics-in-the-cloud
command-line access through the cloud sdk toolkit gsutil: https://cloud.google.com/storage/docs/gsutil

Disclosure: I work at the Broad Institute, which uses a lot of GCP infrastructure; I do not get any compensation from Google Cloud, but I recently wrote a book about Genomics in the Cloud that uses GCP.

I have recently open sourced a Java NIO SPI for S3 that allows direct reads from S3 using Java applications that use the NIO model. This allows applications like GATK and some parts of Picard to understand s3:// URIs.
It is used as a drop in so as long as the jar is on your classpath or in your JRE ext folder then Java will delegate nio operations to this lib when it "sees" an s3:// URI. It will use the credentials and configuration of your current profile so no additional config is require and no recompilation of an application is needed.
The github repo is https://github.com/awslabs/aws-java-nio-spi-for-s3
Disclosure: I am the author of the library, employed by the AWS Health AI Genomics Team, and have financial interests in Amazon


text:
what is the most complete vcf file of population allele frequencies that can be built/downloaded from public datasets?
What is the most complete vcf file of population allele frequencies that can be built/downloaded from public datasets nowadays?
About 5 years or so ago, it used to be the latest release of the CSHL HapMap 12 populations that were part of the 1000 genomes project. These were:
CHB CHD MEX GIH TSI LWK ASW MKK HCB JPT CEU YRI
For example, the EnsEMBL project currently has these 12 populations available as a vcf here:
http://ftp.ensembl.org/pub/current_variation/vcf/homo_sapiens/

CSHL-HAPMAP-HAPMAP-CHB.vcf.gz
CSHL-HAPMAP-HAPMAP-CHD.vcf.gz
CSHL-HAPMAP-HAPMAP-MEX.vcf.gz
CSHL-HAPMAP-HAPMAP-GIH.vcf.gz
CSHL-HAPMAP-HAPMAP-TSI.vcf.gz
CSHL-HAPMAP-HAPMAP-LWK.vcf.gz
CSHL-HAPMAP-HAPMAP-ASW.vcf.gz
CSHL-HAPMAP-HAPMAP-MKK.vcf.gz
CSHL-HAPMAP-HapMap-HCB.vcf.gz
CSHL-HAPMAP-HapMap-JPT.vcf.gz
CSHL-HAPMAP-HapMap-CEU.vcf.gz
CSHL-HAPMAP-HapMap-YRI.vcf.gz

Each SNP has an AF entry from which a multi-populations vcf with rsids, alleles and frequencies can be built, where the AFs are such as AF_CHB, AF_CHD, etc.
The 1000 genomes project populations documentation describes more than these 12 populations, but I haven't seen equivalent population vcfs with AFs built from the individuals within the population for the remainder of these, apart from the original HapMap 12 marked with a * below:
###Populations and codes

*         CHB   Han Chinese             Han Chinese in Beijing, China
*         JPT   Japanese                Japanese in Tokyo, Japan
         CHS    Southern Han Chinese    Han Chinese South
         CDX    Dai Chinese             Chinese Dai in Xishuangbanna, China
         KHV    Kinh Vietnamese         Kinh in Ho Chi Minh City, Vietnam
*         CHD   Denver Chinese          Chinese in Denver, Colorado (pilot 3 only)

         CEU    CEPH                    Utah residents (CEPH) with Northern and Western European ancestry
*         TSI   Tuscan                  Toscani in Italia
         GBR    British                 British in England and Scotland
         FIN    Finnish                 Finnish in Finland
         IBS    Spanish                 Iberian populations in Spain

*         YRI   Yoruba                  Yoruba in Ibadan, Nigeria
*        LWK    Luhya                   Luhya in Webuye, Kenya
         GWD    Gambian                 Gambian in Western Division, The Gambia
         MSL    Mende                   Mende in Sierra Leone
         ESN    Esan                    Esan in Nigeria

*         ASW   African-American SW     African Ancestry in Southwest US
         ACB    African-Caribbean       African Caribbean in Barbados
         MXL    Mexican-American        Mexican Ancestry in Los Angeles, California
         PUR    Puerto Rican            Puerto Rican in Puerto Rico
         CLM    Colombian               Colombian in Medellin, Colombia
         PEL    Peruvian                Peruvian in Lima, Peru

*         GIH   Gujarati                Gujarati Indian in Houston, TX
         PJL    Punjabi                 Punjabi in Lahore, Pakistan
         BEB    Bengali                 Bengali in Bangladesh
         STU    Sri Lankan              Sri Lankan Tamil in the UK
         ITU    Indian                  Indian Telugu in the UK

I presume there is more public data nowadays that can be accessed to build a more complete vcf files of allele frequencies per population, with as many populations as possible to use as a reference dataset with new data.
Thanks in advance.

The best set I'm aware of can be downloaded from https://gnomad.broadinstitute.org .

Thanks, I have seen there is a Downloads section which hopefully will have the files that I am after.

I could only find continent-level AFs for gnomAD. I presume nobody built the population-level frequencies file for these individuals, but maybe someone did and this is available on a website?

Can you explain what you mean by "population-level?" I see ones for racial and other ethnic groups as well.


text:
fixing chromosome with new contig
I want to fix a chromosome of around 50-100Mb (chr.fasta) for a region that is missing a few kilobases of data. I have the missing sequence with a window upstream and downstream that matches perfectly to the reference (contig.fasta).
What's the best way to produce a new chromosome, fixed.fasta that contains the new data from contig.fasta?
See depiction below:
chr.fasta    |---------------------------------------------------   -------------|
contig.fasta                                                   |-===-|

Final intended result:
fixed.fasta  |---------------------------------------------------===-------------|

Any ideas?

I have not received any answer to the question, so I created a small tool to fix the file myself.
I am assuming as per the question that there are perfect matches upstream and downstream to the new sequence, and that anything in between can be deleted/updated by the new sequence.
Given the first_n and last_n bases of the diagram in the question, and the entire patch sequence, we can create a single-line string of the reference with seqtk seq -l9999999 - and then edit it with the line below (in Perl regex):
$chromosome_str =~ s/$first_n.+$last_n/$patch/i

This creates the new string that can be used as the new fasta file.

You can give a try with megamerger from the emboss suite if you have (almost) perfect sequence match.
In case the gap within your chr.fasta complicates the sequence alignment and hence the resulting contig, you can generate a "tiling path" with: upstream of the gap from chr.fasta + contig.fasta + downstream of the gap from chr.fasta (there should be sequence overlaps of course)
In a (more realistic) scenario where there is no perfect sequence match between the overlapping regions, you can fall back to emboss merger, however, it is not designed to work on large sequences.


text:
where to find pathogenicity metadata for E.coli genomes
I have a list of about 900 E.coli genome ids (Genbank plus NCTC ids), e.g.:
NCTC8196_1
NCTC11105_1
CP027464.1
CP027459.1
CP027363.1
CP027388.1
CP027355.1
CP027675.1
CP027445.1
CP027325.1
CP027219.1

and I would like to find pathogenicity metadata for them so that I can classify the molecular data I found with the metadata.
Is there a way to get a table of metadata from these genomes that I can use to associate with my molecular data?
NOTES:
After some searching online, I found Enterobase (http://enterobase.warwick.ac.uk), but it seems the site is mostly dedicated to MLST typing, thus they use different naming conventions than the Genbank id. Maybe there is a way to link the info in Enterobase with Genbank ids somehow?
It seems like Enterobase contains project ids, e.g. PRJNA218110. These project ids can also be found on the other side in the GenBank entries, for which we could associate the GenBank chromosomal ids.
Searching for PRJNA218110 in Enterobase gives a very long list of entries that seem to be from Illumina sequencing data, but in this list, the 'Path/Nonpath' column is always empty.
Enterobase seems to contain lots of entries with useful information, such as: the strain is pathogenic, and produced diarrhea in humans. So a table with entries like the one below would be very useful:
Genbank_id: CPXXXXXXX.1
Path/Nonpath: Pathogenic
Disease: Diarrhea

Genbank_id: CPYYYYYYY.1
Path/Nonpath: Nonpathogenic
Disease: NA

What kind of data, exactly? CP027459.1 is an entire genome. Are you trying to see if this strain is a pathogen? If so, for what species? Are you instead trying to see if mutations in this genome can be pathogenic for the bacterium? Please [edit] and clarify what you're after.

I added a comment now referring to what I can find in Enterobase, thanks @terdon .

Thanks! So you are looking for information on which of these bacteria are pathogens in humans? Humans only or do you also want pathogenicity information on any species?

Humans mainly, but any would be great

you can use BacDive. This is the largest database for bacterial metadata.:
https://bacdive.dsmz.de/
e.g. you can search for the sequence accession "NCTC8196" and you would get the metadata of the according strain:
https://bacdive.dsmz.de/search?search=NCTC8196
What kind of data is available, is highly dependent on the strain. Usually type strains are more rich on data that others.
For further questions don't hesitate to contact us:
contact@bacdive.de


text:
alternatives to PHASTER phage annotation server
I would like to use PHASTER (http://phaster.ca) from the University of Alberta, but it's currently down for maintenance, as is its predecessor, PHAST.
What alternatives to PHASTER are there for searching and annotating phage sequences within assembled genomes?
EDIT:
I have seen that PROKKA has some entries mentioning the possibility of using PHASER (currently discarded as an option) but also people using their own HMMs to annotate phage sequences using PROKKA, e.g.:
https://github.com/tseemann/prokka/issues/375
Any ideas?

Phigaro https://github.com/bobeobibo/phigaro, Virsorter

Thanks @Chris_Rands , I had a look at phigaro, and I like the fact that the pVOG families give an annotation for things like tail fiber proteins and the rest. Useful.

The answer is to use phastaf, a recently developed software tool developed by Torsten Seemann (the same developer as prokka), which comes with a packaged phaster.ca db.  Get it here

Thanks, I will have a look.


text:
IGV jump to rsID location
Is it possible to load a list of rsID identifiers (e.g. rs66823261) to IGV to then be able to jump to the location of each of them when copy+pasting the rsID on the locus text box?
For human and other genomes, there is currently the feature of typing in a gene name, then jumping to the gene location.
I wonder if there is a way of achieving the same for rsIDs.
EDIT: thanks for the answers so far. I downloaded the Infinium GSA bed file from Illumina's website (https://support.illumina.com/downloads/infinium-global-screening-array-v1-0-support-files.html), loaded on an IGV session, which asked to create an index. I did a liftover of the file to hg38 and sorted with bedtools sort. Then I went to the text box, introduced one of the rs ids:
rs6569648
Got a positive result. Thank you all!

Yes, if you load the list you are interested in as a non-indexed "BED" file then you can jump to an rsID or any other name in the BED file. You'll need a targeted list, dbSNP for example is much too large to load in this way.

Great, thank. Is there an hg38 bed file with the rsIDs of one of the 1M SNP arrays out in the market?

you can use mysql ucsc to get the position of the SNP and then ask IGV to take some picures: see https://software.broadinstitute.org/software/igv/batch
$ echo -e 'rs66823261\nrs25' | awk '{printf("select concat(chrom,\":\",chromStart-10,\"-\",chromEnd+10) from  snp150 where name=\"%s\";\n",$1);}' | mysql --user=genome --host=genome-mysql.soe.ucsc.edu -A -D hg38 -N | awk '{printf("goto %s\nsnapshot\n",$1);}'

goto chr8:170681-170702
snapshot
goto chr7:11584131-11584152
snapshot

If you want all rsIDs (complete dbSnp database) then Pierre's answer is the best option currently.    Any bed file with all the rsIDs would be too big to load.  For a targeted list of, say, 100,000 or fewer snps you can use a bed file.
BTW igv.js already does the equivalent of Pierre's suggestion, this will be ported to IGV desktop in the near future.


text:
De novo antibody sequencing fromMS/MS Ion Trap proteomics raw signal
I've read about different proteomics tools that can do de novo assembly of MS/MS Ion Trap proteomics signatures for de novo antibody sequencing. Given the way CDR3s of antibodies are created, together with Somatic Hypermutation, can generate a huge diversity in the CDR regions of the heavy and light chains, the reference-based methods can fall short of determining the full sequence of an unknown antibody without a set of reference sequences including those with its CDRs.
Are there any other tools that can determine antibody sequences from MS/MS proteomics raw signals?

A few notes on data generation first:
There are currently two broad approaches for the generation of bottom-up or "shotgun" MS proteomics: data-dependent acquisition (DDA) and data-independent acquisition (DIA).1 In tandem MS (MS/MS), the DDA approach only puts forward certain peptides generated during the first cycle of MS for fragmentation during the second cycle, while with the DIA approach, all peptides generated during the first MS cycle can be fragmented in the second round.
The recently-introduced Orbitrap Fusion mass spectrometry permits various types of MS2 acquisition methods. There are four main MS2 strategies:

HCD-OT (higher-energy-collisional-dissociation with Orbitrap detection),

HCD-IT (HCD with ion trap, IT),

CID-IT (collision-induced-dissociation with IT) and

CID-OT on Orbitrap Fusion

Reading from the biorxiv preprint of Stitch:
https://www.biorxiv.org/content/10.1101/2022.03.07.483237v1
"With the aid of specialized software packages like DiPS, Supernovo, or ALPS, full and accurate sequences of the heavy and light chains can be reconstructed with the MS/MS spectra of the digested peptides."
Before running Stitch or similar tools, a peptide profiling tool needs to be executed on the .raw or .mgf data alone.
STITCH: https://www.biorxiv.org/content/10.1101/2022.03.07.483237v2
This currently requires the output of analysing .raw data via PEAKS, which is not free unfortunately. Authors currently looking into using other free de novo sequencing software for input into Stitch. The software also allows for FASTA format, which can then be executed on the output of PepNovo, Novor or similar.
DiPS: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5461544/
ALPS: https://www.nature.com/articles/srep31730
"The results show that ALPS was able to assemble three complete monoclonal antibody sequences of length 216–441 AA, at 100% coverage and 96.64–100% accuracy."
DeepNovo: https://www.pnas.org/doi/10.1073/pnas.1705691114
"This field has been actively studied over the past 20 y, and many de novo sequencing tools have been proposed, such as PepNovo, PEAKS, NovoHMM, MSNovo, pNovo, UniNovo, and Novor among others (6–19). The recent "gold rush" into mAbs has undoubtedly elevated the application of de novo sequencing to a new horizon (20–23)"
"DeepNovo considerably outperformed state of the art methods, achieving 7.7–22.9% higher accuracy at the amino acid level and 38.1–64.0% higher accuracy at the peptide level"
Available for download, requires activation via sending email to author, which then points to the company behind PEAKS, which will give you a 15 day evaluation license.
NovoHMM: the download page leads to a missing URL
DeepResolution2: https://www.sciencedirect.com/science/article/abs/pii/S0039914022002119
Capable of analyzing GC-MS data rather than MS/MS Ion Trap
DePS: https://arxiv.org/abs/2203.08820
DeepNovoV2: https://github.com/volpato30/DeepNovoV2
DPST: https://arxiv.org/abs/2203.13132
Link does not work https://github.com/Yan98/DPST
PepNovo: http://proteomics.ucsd.edu/software-tools/531-2/
The new version is from 2012, includes models for CID and LTQ, but not HCD models.
pNovo: http://pfind.org/
Includes models for CID, HCD, ETD and HCD+ETD.
Novor: https://www.rapidnovor.com/
Filling in a form, someone will contact you. Doesn't seem to be free.
METHODOLOGY used so far:
1 - ThermoRawFileParser: Convert .raw data to .gcf
2 - PepNovo: produce list of peptides, one-liner to convert to .fasta
3 - Stitch: analyse given a list of V(D)J+C references from templates and the .fasta file from above.