text:
samtools mpileup empty when filtering out flags
I produced a bam file by aligning reads to a small set of synthetic sequences using bwa-mem.
I am heavily filtering reads that are not paired and of a certain orientation.
Applying the filtering, I get a few thousands of reads:
samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools view -h -f32 -  | wc -l

Gives me 89502 reads.
If I then pipe this into samtools mpileup, I get no results:
samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools view -h -f32 -  | \
samtools mpileup --excl-flags 0 -Q0 -B -d 999999 - | wc -l

Returns 0.
I tried different combinations of filtering, and when I do both -f 16 and -f 32 returns empty, but if I do either of those, then it works:
samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools mpileup --excl-flags 0 -Q0 -B -d 999999 - | wc -l

Returns 1056.
Any ideas why? My thinking was that it would work with --excl-flags 0.
EDIT: substituting mpileup for depth does work, and prints out each position and the depth as expected.
EDIT2: adding -q 0 to mpileup gives the same empty result.
Thanks in advance

Hi, could you add some details of how you got the synthetic reads and how you have aligned the to the reference?

Are you aware that -f 16 combined with -f 32 (btw, you can just -f 48 and skip the piping) are filtering for discordant reads? What happens if you add -q 0 to mpileup?

Since depth and mpileup share a fair bit of code, I have to wonder if you simply have no variants there. What happens if you specify -a and filter out the positions with no coverage?

BTW, I suspect this will be a case where you'll need to post the BAM file somewhere. There are a few (currently commented out) printf()s in the samtools code that could probably be used to track down what's actually happening here.

Devon already mentioned it in passing but it’s too important to get lost: Samtools flags are bit masks that you need to combine! Don’t call samtools repeatedly; call it once, with the combined bit flag. You can use the BROAD online tool to compute the combined flag if you don’t want to do it manually.

By using -h in the samtools view command, you're including all the header lines in your word count. If you happen to have about 89500 reference sequences, then the lengths of those would all appear in the header and inflate the -h word count, but not the mpileup count. Try piping it through an additional samtools view (i.e. without -h) and see if the counts change:
...
samtools view -h -f32 -  | \
samtools view | wc -l

Also, samtools mpileup by default only considers high-quality bases and concordant reads. Try adding a -A to your mpileup line (which stops anomalous read pairs from being discarded):
...
samtools mpileup -A -Q0 -B -d 999999 - | wc -l

Whether or not this is actually a good idea will be dependent on what you want to get out of the analysis, and what the downstream programs / analyses are expecting.

Bingo!  Adding -A, --count-orphans     do not discard anomalous read pairs did the trick!

i do have a follow up question. Is using the -A option a good idea? I have same problem as OP (i.e. empty mpileup), but can't find any useful info. of course with -A it seems to work but...how reliable will all downstream analysis be?

As with most questions of this nature, there's not a universal answer. If you want a useful answer, please post this as a separate question with enough context to explain what the eventual end-goal is, and what that downstream analysis is.


text:
how to download all metadata for a single TCGA dataset and link to data file UUIDs?
I've been able to download specific data types from the TCGA portal (https://portal.gdc.cancer.gov/) using the Manifest files and the GDC Tool.
But it's not clear to me how to download the metadata (e.g. cancer type, stage, etc.) with the links to the UUIDs of a specific data type, e.g. the UUID of a methylation data file to the case UUID:
https://portal.gdc.cancer.gov/files/14dc6f1f-21d0-4b05-9ad5-1d52f5020b12
Any ideas?

To download metadata associated with a specific data file, such as the cancer type and stage, you can use the GDC API (Application Programming Interface). The API allows you to query the GDC database and retrieve information in a structured format, such as JSON or XML.
To download the metadata for a specific data file, you will need the file's unique identifier, also known as its UUID. In your example, the UUID of the methylation data file is "14dc6f1f-21d0-4b05-9ad5-1d52f5020b12".
To use the GDC API, you can construct a query using the GDC API endpoint and the file's UUID. For example, to retrieve the metadata for the methylation data file with the UUID "14dc6f1f-21d0-4b05-9ad5-1d52f5020b12", you could use the following query:
https://api.gdc.cancer.gov/files/14dc6f1f-21d0-4b05-9ad5-1d52f5020b12?expand=metadata
This query uses the "files" endpoint of the GDC API and specifies the file's UUID. The "expand=metadata" parameter tells the API to return the metadata associated with the file. The response will be in json format:
{"data": {"data_format": "TXT", "access": "open", "file_name": "5b6353ee-af39-4bb4-a801-1454dfbc338b.methylation_array.sesame.level3betas.txt", "submitter_id": "d67719e6-4149-4492-a265-1f5d5689b904", "data_category": "DNA Methylation", "acl": ["open"], "type": "methylation_beta_value", "platform": "Illumina Human Methylation 450", "file_size": 13168876, "created_datetime": "2022-02-15T17:17:23.357435-06:00", "md5sum": "ab0c90db97a3e12634ef3f54f3ec835d", "updated_datetime": "2022-02-16T19:13:44.087132-06:00", "file_id": "14dc6f1f-21d0-4b05-9ad5-1d52f5020b12", "data_type": "Methylation Beta Value", "state": "released", "experimental_strategy": "Methylation Array", "version": "1", "data_release": "32.0 - 36.0"}, "warnings": {"expand": "unrecognized values: [metadata]"}}

source: Stack Exchange

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

sort: votes


text:
low-memory high-speed replacement for Picard MarkDuplicates
I am running Picard MarkDuplicates with the following parameters below. On the file described, it takes about 41.6Gb of RAM memory and about 20-25 minutes to compute (only uses 1 core AFAICS).
java -Xmx54890m -jar picard.jar MarkDuplicates \
    I=/home/dnanexus/in/sorted_bam/TST73-86-2IC_S16_L00.bml.GRCh38.karyo.bam \
    O=./TST73-86-2IC_S16_L00.bml.GRCh38.karyo.deduplicated.bam \
    M=./TST73-86-2IC_S16_L00.bml.GRCh38.karyo.duplication_metrics \
    CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT REMOVE_DUPLICATES=true

What equivalent free/open source software (samtools?, sambamba?) can I use to obtain the same deduplicated bam file but at a lower cost, considering the instance used here has 8 cores but only one of them is being used, and has about 68Gb of RAM out of which it's using 46Gb odd at peak.
I think previous versions of samtools/sambamba didn't behave the same way as Picard MarkDuplicates with the parameters above, but I would like to know which recent versions would give me the same output.

I assume you only want free (in terms of price, at least) software, right? Sentieon have a tool for this that is orders of magnitude faster (and I would guess uses less RAM but I'd have to check) but you need to pay for it.

Correct, open source.

I done no benchmark, but I guess samblaster is bit faster (by my subjective perception of time).

You should specify a smaller "-Xmx". MarkDuplicate doesn't need that much memory, but seeing a large "-Xmx", java will use extra.

As you suggested, sambamba is faster at marking duplicates than picard (it's also multithreaded). Recent versions of samtools have a rewritten duplicate marking algorithm, though I doubt it'll be as quick as sambamba. 46GB of RAM seems excessive for marking duplicates unless you're having it store the whole file in memory.

Does sambamba with default parameters equate to Picard MarkDuplicates?

I'm not sure that any tool exactly equates (there are a lot of edge cases in duplicate marking and picard misses some things), but it's a sufficient replacement. Sambamba is also used for this purpose by a number of international consortia (I think either DEEP or IHEC, if I remember correctly).

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


text:
What do the FASTQ file names mean here?
I would like to run an analysis on Nanopore data. I'm trying to download sample data from https://github.com/nanopore-wgs-consortium/chm13. I would like to run minimap2 on the FASTQ files, resulting in a single BAM file.

Q: What do the file names mean? What is "S13"? What is "L002_l1"? What is "L002_R1"? Am I supposed to download them and combine them into a single FASTQ file, and then align them?

What is "S13"?
S# is most likely the sample number. If so, you do not want to combine these.
What is "L002_l1"?
L00# is the lane number (of the 8 Illumina flow cell lanes). 'I1_' (not 'L') are index files and are not needed (see also https://www.biostars.org/p/374581/).
What is "L002_R1"?
R1 is one end of paired-end sequencing, R2 is the other end. These filepairs must have the same number of lines. Use an aligner that utilizes both.

Peter's answer is a good one and to the point. I'll be a bit more verbose.
The breakdown of FASTQ file names that come directly from the sequencer typically have the following format:
{sample_name}_S{sample number}_L{lane number}_{R/I}{read or index number}_001.fastq.gz
sample_name is the sample name provided by you (or whoever sequenced the data) to the sequencer.
To save time and money, samples are often "multiplexed" on sequencing runs. This multiplexing means you aggregate samples ("S") on the same run and split samples on different lanes ("L") to limit some technical noise. Lanes are essentially technical replicates for the sequencing data on the same sequencing run, and these numbers tell you which sample and which lane that FASTQ comes from.
Sequencing can be done in single-end or paired-end. Single-end sequencing generates a single FASTQ file (_R1), whereas paired-end sequencing produces 2 FASTQ files (_R1 and _R2).
The data you highlight isn't actually Nanopore:

Approximately 50x of data was generated on a NovaSeq instrument.

10x includes "indexes" ("I") in addition to the sample and lane information to keep track of things like single cells or other unique molecular identifiers, and that information is listed after the I. _R1, _R2, and _I1 all go together for a unique sample and lane combination

The file names mean this isn't Nanopore data. Read the bit just above these file names:

Approximately 50x of data was generated on a NovaSeq instrument.
Based on the summary output of Supernova, there are 1.2 billion reads with 41x effective coverage.

That linked github page has Nanopore, Illumina/10x, BioNano, HiC, and PacBio data. For the nanopore files, you should be looking at this section.

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


text:
How to use Python to count k-mers?
I have some FASTQ sequence files and a FASTA file for some regions I'm interested in.
I would like:

Build an index for the FASTA file
Use the index to count number of k-mers occurred in my sequence files

I know how to do this in many k-mer counting tools, such as Kallisto, Jellyfish, DSK, Salmon. However, I'm more interested in Python (khmer is a Python package, all the other k-mer counting tools are C++).
Unfortunately, I can't figure out how to do that in khmer from the documentation:

http://khmer.readthedocs.io/en/v2.1.1/introduction.html

It talks about RNA-Seq, metagenomics, de Bruijn graph etc... All I need are just counts. Khmer is a popular k-mer tool, so it should be able to do it, but how?
Q: How to build an index and count k-mers using khmer?
EDIT:
I don't need number of unique k-mers, I'm looking for a solution like Jellyfish where we build an index and count k-mers to the index.
Something like:
MY_SEQ_1  100 counts
MY_SEQ_2  134 counts

Does this accomplish what you're looking for: https://github.com/dib-lab/khmer/blob/master/examples/python-api/exact-counting.py

If you want to count the number of unique k-mers that occur in your data set, you should use the unique-kmers.py script, which implements a HyperLogLog-based cardinality estimator.
If the number of unique k-mers is not what you're after, please clarify.
EDIT:
khmer uses probabilistic data structures internally, which store k-mers as hashed values that cannot be un-hashed back into k-mer sequences uniquely. Querying a k-mer's abundance requires you to know which k-mer(s) you're looking for, which in this case would require a second pass over the reads.
The closest thing to what you're asking provided in khmer's command-line scripts is abundance-dist.py (or the alternative abundance-dist-single.py), which will produce a k-mer abundance histogram but not per-kmer abundance.
The load-into-counting.py scripts computes the k-mer abundances and stores them in a probabilistic data structure, which is written to disk and can subsequently be re-loaded into memory quickly. Many of the scripts in khmer expect the reads to have been pre-processed by load-into-counting.py, while others will invoke the counting routines directly themselve. It's also fairly easy to do this with the Python API. In every case, the accuracy of the k-mer counts is dependant on memory considerations.
>>> import khmer
>>> counts = khmer.Counttable(31, 1e7, 4)
>>> counts.consume_seqfile('reads.fq.gz')
(25000, 1150000)
>>> counts.get('TGACTTTCTTCGCTTCCTGACGGCTTATGCC')
117

To check the false positive rate:
>>> fpr = khmer.calc_expected_collisions(counts)
>>> fpr

If you must have the counts of each k-mer, it's not much more work...although reporting each k-mer's abundance only once would consume a lot of memory with a naive approach.
>>> outfile = open('outfile.txt', 'w')
>>> seenkmers = set()  # Consumes a lot of memory for large input!!!
>>> for read in khmer.ReadParser('reads.fq.gz):
...     for kmer in counts.get_kmers(read.sequence):
...         if kmer not in seenkmers:
...             print(kmer, counttable.get(kmer), file=outfile)
...             seenkmers.add(kmer)

Not exactly what I'm looking for. EDITED.

See this question if you run into problems not counting kmers above 255

I wrote a command-line k-mer counter called kmer-counter that will output results in a form that your Python script can consume: https://github.com/alexpreynolds/kmer-counter
You can grab, build and install it like so:
$ git clone https://github.com/alexpreynolds/kmer-counter.git
$ cd kmer-counter
$ make
$ cp kmer-counter /usr/local/bin

Once the binary is in your path, you might use it in Python like so:
k = 6
fastaFile = '/path/to/some/seqs.fa'
kmerCmd = 'kmer-counter --fasta --k=%d %s' % (k, fastaFile)
try:
    output = subprocess.check_output(kmerCmd, shell=True)
    result = {}
    for line in output.splitlines():
        (header, counts) = line.strip().split('\t')
        header = header[1:]
        kmers = dict((k,int(v)) for (k,v) in [d.split(':') for d in counts.split(' ')])
        result[header] = kmers
    sys.stdout.write("%s" % (str(result)))
except subprocess.CalledProcessError as error:
    sys.stderr.write("%s" % (str(error)))

Given example FASTA like this:
>foo
TTAACG
>bar
GTGGAAGTTCTTAGGGCATGGCAAAGAGTCAGAATTTGAC

For k=6, you would get an iterable Python dictionary like this:
{'foo': {'TTAACG': 1, 'CGTTAA': 1}, 'bar': {'GTTCTT': 1, 'AGAACT': 1, 'GAGTCA': 1, 'ATGGCA': 1, 'GAACTT': 1, 'ATTCTG': 1, 'CTAAGA': 1, 'CTTCCA': 1, 'ATTTGA': 1, 'GGAAGT': 1, 'AGGGCA': 1, 'CCTAAG': 1, 'CTCTTT': 1, 'AATTTG': 1, 'TCTGAC': 1, 'TTTGCC': 1, 'CTTAGG': 1, 'TTTGAC': 1, 'GAAGTT': 1, 'CCCTAA': 1, 'AGAATT': 1, 'AGTCAG': 1, 'CTGACT': 1, 'TCTTAG': 1, 'CGTTAA': 1, 'GTGGAA': 1, 'TGCCAT': 1, 'ACTCTT': 1, 'GGGCAT': 1, 'TTAGGG': 1, 'CTTTGC': 1, 'TGGAAG': 1, 'GACTCT': 1, 'CATGCC': 1, 'GCAAAG': 1, 'AAATTC': 1, 'GTCAAA': 1, 'TGACTC': 1, 'TAGGGC': 1, 'AAGTTC': 1, 'ATGCCC': 1, 'TCAAAT': 1, 'CAAAGA': 1, 'AACTTC': 1, 'GTCAGA': 1, 'CAAATT': 1, 'TAAGAA': 1, 'CATGGC': 1, 'AAGAAC': 1, 'AAGAGT': 1, 'TCTTTG': 1, 'TTCCAC': 1, 'TGGCAA': 1, 'GGCAAA': 1, 'AGTTCT': 1, 'AGAGTC': 1, 'TCAGAA': 1, 'GAATTT': 1, 'AAAGAG': 1, 'TGCCCT': 1, 'CCATGC': 1, 'GGCATG': 1, 'TTGCCA': 1, 'CAGAAT': 1, 'AATTCT': 1, 'GCATGG': 1, 'ACTTCC': 1, 'TTCTTA': 1, 'GCCATG': 1, 'GCCCTA': 1, 'TTCTGA': 1}}

You can use standard Python calls to manipulate this dictionary object and get sums of counts per record, for sequence, etc. which seems to answer your question. Please feel free to clarify what you're looking for if this object representation is not clear.
For a fully-fleshed out demonstration, see: https://github.com/alexpreynolds/kmer-counter/blob/master/test/kmer-test.py

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


text:
The effects of incomplete bisulfite conversion upon mapping efficiency
This question has also been posted on Biostars
I have sequenced numerous multiplexed pools of BS amplicon-seq libraries derived from human samples on a MiSeq over the past few weeks. I have been utilising trim-galore and Bismark for alignment and am finding the mapping efficiency to be really low for two pools (55% & 30% respectively) both of which had a cytosine per base sequence of around 10-20% throughout the entire read in their fastqc files.
The high C content visible in the fastqc files makes me think the poor mapping efficiency is due to poor bisulfite conversion, as this would be expected to be close to zero is bisulfite conversion had actually took place.
I am new to the field so any help would be greatly appreciated.

It may be good to familiarize with the steps behind mapping bisulfite converted reads. In this post on biostars I give a quick breakdown of the steps involved if you want to have a look (but of course there is a lot of material available on it).

Bisulfite conversion efficiency has no effect on the mapping rate in bismark and similar tools. The reason is that the reads are fully bisulfite converted in silico before alignment to minimize mapping bias.
I would suggest that you play around with the settings handed to bowtie2, such as using local alignment and modifying the --score-min option to allow more mismatches. Alternatively, you might try a different aligner like bwameth, which will always do local alignment.

Thank you. --local flag in the bowtie2 docs looks promising, I'll give this a shot when I get into the office. Will Bismark parse any bowtie2 options/flags to bowtie2 (I ask as --local is not present in bisamark docs but only within bowtie2 docs)? What can I conclude about the high cytosine per base sequence in the fastqc files? I've generally found this is at 0% in previous analyses or 10% at the ends, at which point i would trim the ends.

I don't think that bismark will handle the --local flag correctly, since it expects reads to either fully align or be discarded, without any soft clipping.

If after running bismark the MBIAS plots don't look relatively flat, you can pass a filter parameter to bismark_methylation_extractor, so that e.g. you ignore the last 10bp of each read.

@719016: You might be correct, it's been a few years since I've looked at the bismark code, I generally try to avoid using it since it's so incredibly slow and used to give poor results.

As Devon Ryan also proposes, you may want to try bwameth, which is a simple wrapper around bwa-mem, and does keep the soft-clipped alignments if they exist. You can then apply another script for counting the CpGs from the bam, e.g. this one: https://github.com/dariober/bioinformatics-cafe/tree/master/bam2methylation

I generally recommend MethylDackel for methylation extraction and QC, but I'm biased :P


text:
Obtaining uniquely mapped reads from BWA mem alignment
This is based on a question from betsy.s.collins on BioStars. The original post can be found here.

Does anyone have any suggestions for other tags or filtering steps on BWA-generated BAM files that can be used so reads only map to one location? One example application would be to find seeds for the TULIP assembler/scaffolder, which works best for reads that map to unique genomic locations.

The questioner is referring to genomically-unique alignments, which are something different from multiple possible alignments at a single location. If the only multiple alignment is at the same genomic location, those alignments are still unique in the genomic sense. The impossibility of finding the exact, correct, alignment is a well-known problem, and there are a few downstream methods (e.g. left normalisation) to make sure that multiple local alignments and/or sequencing errors have reduced impact.
The concept of "uniquely mapped reads" is a loaded term, and most sources suggest filtering by MAPQ should do the trick. However, this approach doesn't seem to work when using BWA as a read mapper.
Uniquely mapped reads (i.e. a read that maps to a single location in the genome) are sometimes preferred when running analyses that depend on the quantification of reads, rather than just the coverage (e.g. RNASeq). Duplicated reads require additional filtering or normalisation on the analysis side, and most downstream programs won't consider concepts like "a fifth of the read maps somewhere here on chromosome 5, and four-fifths of the read maps somewhere here on chromosome 14."

Update - as of January 2021, samtools can now do filtering based on an expression that includes tag variables. In this case, this expression can be used to exclude any reads that have either an XA or SA tag:
samtools view -b mapped.bam -e '!([XA] | [SA])' > unique_mapped.bam

For more details on the samtools expression parameter, see the samtools documentation:
https://www.htslib.org/doc/samtools.html#FILTER_EXPRESSIONS

Original answer follows....
To exclude all possible multi-mapped reads from a BWA-mapped BAM file, it looks like you need to use grep on the uncompressed SAM fields:
samtools view -h mapped.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > unique_mapped.bam

Explanation follows...
I'm going to assume a situation in which a bioinformatician is presented with a mapped BAM file produced by BWA, and has no way of getting the original reads. One high-effort solution would be to extract the mapped reads from the BAM file and re-map with a different mapper that uses the MAPQ score to indicate multiple mappings.
... but what if that were not possible?
My understanding of BWA's output is that if a read maps perfectly to multiple genomic locations, it will be given a high mapping quality (MAPQ) score for both locations. Many people expect that a read that maps to at least two locations can have (at best) a 50% probability of mapping to one of those locations (i.e. MAPQ = 3). Because BWA doesn't do this, it makes it difficult to filter out multiply-mapped reads from BWA results using the MAPQ filter that works for other aligners; this is likely to be why the current answer on Biostars [samtools view -bq 1] probably won't work.
Here is an example line from a BWA mem alignment that I've just made. These are Illumina reads mapped to a parasite genome that has a lot of repetitive sequence:
ERR063640.7     16      tig00019544     79974   21      21M2I56M1I20M   *       0       0       TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC    <AABBEGABFJKKKIM7GHKKJK>JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG>D    NM:i:4  MD:Z:83A13      AS:i:77 XS:i:67 XA:Z:tig00019544,-78808,21M2I56M1I20M,6;tig00019544,-84624,79M1I20M,6;tig00019544,-79312,33M4I42M1I20M,8;

BWA mem has found that this particular read, ERR063460.7, maps to at least three different locations: tig00019544, tig00019544, and tig00019544. Note that the MAPQ for this read is 21, so even though the read maps to multiple locations, MAPQ can't be used to determine that.
However, the alternative locations are shown by the presence of the XA tag in the custom fields section of the SAM output. Perhaps just filtering on lines that contain the XA tag will be able to exclude multiply-mapped reads. The samtools view man page suggests that -x will filter out a particular tag:
$ samtools view -x XA output.bam | grep '^ERR063640\.7[[:space:]]'
ERR063640.7 16  tig00019544 79974   21  21M2I56M1I20M   *   0   0   TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC    <AABBEGABFJKKKIM7GHKKJK>JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG>D    NM:i:4  MD:Z:83A13  AS:i:77 XS:i:67

... so it filtered out the tag (i.e. the tag no longer exists in the SAM output), but not the read. There are no useful bits in the FLAG field to indicate multiple genomic mappings (which I know can be filtered to exclude the read as well), so I have to resort to other measures.
In this particular case, I can use grep -v on the uncompressed SAM output to exclude alignment lines that have the XA tag (and re-compress to BAM afterwards, just to be tidy):
$ samtools view -h output.bam | grep -v 'XA:Z:' | samtools view -b > output_filtered.bam
$ samtools view output_filtered.bam | grep '^ERR063640\.7[[:space:]]'
<no output>

Hurray! reads filtered. As a little aside, this grep search has a fairly substantial computational load: it's looking for some string with the text XA:Z: somewhere in the line, and doesn't actually capture every situation. Some masochistic person might come along at a later date and decide that they're going to call all their reads HAXXA:Z:AWESOME!:<readNumber>, in which case a tweak to this grep search would be needed to make sure that there's a space (or more specifically a tab character) prior to the XA:Z:.
Now I do a check for any duplicated read names, just to be sure:
$ samtools view output_filtered.bam | awk '{print $1}' | sort | uniq -d
ERR063640.1194
ERR063640.1429
ERR063640.1761
ERR063640.2336
ERR063640.2825
ERR063640.3458
ERR063640.4421
ERR063640.4474
ERR063640.4888
ERR063640.49
ERR063640.4974
ERR063640.5070
ERR063640.5130
ERR063640.5300
ERR063640.5868
ERR063640.6116
ERR063640.6198
ERR063640.6468
ERR063640.6717
ERR063640.6797
ERR063640.7322
ERR063640.750
ERR063640.7570
ERR063640.7900
ERR063640.8115
ERR063640.8405
ERR063640.911
ERR063640.9206
ERR063640.9765
ERR063640.9986

Oh... damn. I wonder what they are:
$ samtools view output_filtered.bam | grep '^ERR063640.3458[[:space:]]'
ERR063640.3458  16  tig00002961 5402    60  58S38M  *   0   0   AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATATCCGTATCGCGAAGATCAGGACTTACTCCGCAGAAGAA    DD6HFFJBKFH=KDILKLGLJEKLKGFJIH8IKHLLMJEK:L:HBGJIHJKFLLKIHJDHLNKCK;KMKGMFKJILIIIMKI9JLKKHEJFII?CC    NM:i:0  MD:Z:38 AS:i:38 XS:i:0  SA:Z:tig00002377,202353,-,14M3I5M1I35M38S,19,5;
ERR063640.3458  2064    tig00002377 202353  19  14M3I5M1I35M38H *   0   0   AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATA  DD6HFFJBKFH=KDILKLGLJEKLKGFJIH8IKHLLMJEK:L:HBGJIHJKFLLKIHJ  NM:i:5  MD:Z:5G48   AS:i:35 XS:i:27 SA:Z:tig00002961,5402,-,58S38M,60,0;

Aha! Supplemental alignments, which use the official SA tag [other canonical alignments in a chimeric alignment]. These appear to be situations where a single read has been split up, and maps to multiple locations. Note that in this case, both of the alignments still have MAPQ scores of over 3. It sounds like the questioner would also want to get rid of these situations as well. This time, there are standard flag fields as well to deal with these situations (0x800: secondary alignment). Except it's not enough to just filter the supplemental alignment, because both read mappings should be removed, rather than just the one (or ones) that happened to be tagged as secondary.
Luckily, BWA appears to put the SA tag into all reads containing supplementary alignments (if this is not the case, I'm sure someone will correct me on that). So, I add in the SA search as an additional grep filter:
$ samtools view -h output.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > output_filtered.bam
$ samtools view output_filtered.bam | awk '{print $1}' | sort | uniq -d
<no output>

Done. Easy peasy! </sarcasm>
... that original "high-effort" solution of using a different aligner doesn't look so bad now.

that's a neat solution, and well argued. I don't think it's that high effort if you condense the whole thing to a few lines of shell script either. :)

The thing with using a MAPQ to filter out multimappers is that it uses the presumption of a multimapper being "equally well mapping", which for even the examples presented here won't be the case. I've set your reply as the accepted answer on biostars but asked op there to clarify a bit on what exactly is wanted.

Yes, it's not hard if you know what to expect, which is why I've tried to run through all the problems that might happen.

About the read names possibly containing "XA:Z:", a way to achieve robustness would be to use a proper SAM format parser, such as pysam, if one programs in python, rather than grep. I don't know how much slower this would be.

Worst case, you can write a couple lines of python or likely use Pierre's jvarkit to filter things out.

For much faster and stable$ implementation, sambamba can handle this using following one-liner:
sambamba view -t 12 -h -f bam -F "mapping_quality >= 1 and not (unmapped or secondary_alignment) and not ([XA] != null or [SA] != null)" test.bwa-mem.bam -o test-uniq.bam

See manpage and sambamba synatx for filters for details on further usage.
$: Besides grep based search being slow, it may return false positive negative hits (@gringer's reply above). I also do not find awk based approach reliable because fields like XA, SA, etc. are optional fields in SAM format, and not positionally fixed fields.

In what way were the hits false positive?

Hmm.. I guess it should be false negative because you were excluding those search patterns with grep -v, right?

The only likely way grep could produce the wrong results is if someone modified the read names or a read group to contain one of those strings, which seems quite unlikely.

agree, it is quite unlikely but nonetheless, grep based approach is inherently slow.

The slow part is in the conversion from BAM to SAM by samtools view. The grep program is very quick at what it does.


text:
How much does Nanopore cDNA Sequencing Cost?
[this question is based on a question that was asked on Reddit]
We're interested in doing whole-transcriptome cDNA sequencing of 24 mouse cell lines, and are deciding between Illumina and Nanopore sequencing. What are the costs involved in Nanopore sequencing, and how much hassle is it to analyse Nanopore data?

The standard PCR-cDNA barcoding kit is \$750 USD, and has enough reagents for six runs, with each run using up to 24 samples (i.e. 144 samples total). Each run also needs one MinION flow cell (\$900 USD).
The total cost including external reagents is naturally dependent on the number of samples, and gets cheaper as the number of samples increases. For 24 samples, it's about \$5200 NZD (\$3100 USD). See below for more details.
If you could outsource this to a lab that has a PromethION (or P2 Solo), the sequencing would be a great lot cheaper. PromethION flow cells start at \$3600 for a pack of 4 (i.e. the same unit cost as MinION flow cells), but have about ten times the yield of MinION flow cells (which means more samples, or greater depth).
From a bioinformatics perspective, there's more work to do with nanopore reads because there are no commonly used / established pipelines (so I've had to create my own), but the results are much cleaner once they've been processed.
Cost Breakdown
For cDNA sequencing in Aotearoa, MinION works out cheaper than Illumina when compared against a service centre, about two thirds the price for us. Here is the reagent list I wrote up a couple of months ago for sequencing from 30 whole transcriptome samples, mouse or human (prices in NZD). This is with an ideal read threshold of 1M reads per sample (which I have found to give good quantitative results for differential expression). You could probably get away with lower yield for microbial sequencing:
ONT reagents:

1 x PCR-cDNA Barcoding Kit 24 \$750 USD [\$1210 NZD]
1 x Rapid Adapter Auxiliary V14 \$288 USD [\$481 NZD]
1 x R10.4.1 flow cell @ \$900 USD [\$1503 NZD]

Total ONT cost: \$3194
Other reagents:

10 mM dNTP NEB N0447 [\$145 NZD]
Maxima H Minus RT ThermoFisher EP0751 [\$346 NZD]
RNAclean or AmPure XP beads [AmPure XP: \$800 NZD]
Exonuclease 1 NEB M0293 [\$160 NZD]
RNAseOUT ThermoFisher 10777019 [\$569 NZD]

Total external reagent cost (excludes nuclease-free water, tubes, pipette tips, gloves, ethanol, DNA QC): \$2020 NZD
[As well as whatever you need to extract 200ng total RNA per sample]
Our current protocols for this PCR-cDNA barcoding + V14 RA approach can be found here.
The comparative Illumina cost according to the service centre we asked a couple of months ago was \$185 NZD per sample + \$5000 NZD per flow cell (or $9440 NZDfor 24 samples).
For almost all of these external reagents the amount ordered is an absolute minimum (from memory I think most are something like 50 reactions), so doing a second equivalent experiment would mostly only need the flow cells to be purchased.
In addition to the cost, turnaround time with MinION is quicker (for us it's a couple of weeks versus a couple of months for Illumina [with most of that time spent waiting in a queue]), and the results are nicer to analyse: full length transcript sequences with very even coverage profiles.
I'm aware that sequencing is more expensive in Aotearoa than it is in other countries, so the cost / time difference may not be as extreme elsewhere.
Bulk discounts for ONT flow cells
The prices above are for ordering reagents and flow cells one at a time, using the most expensive reagents. I call this a "grant funding" budget, because it's the quote I would give to people who are looking to fund research. There are ways to reduce costs when organising sequencing in bulk, and I've already hinted above at the extra capacity in reagent kits.
However, the biggest savings for nanopore sequencing can be achieved by ordering lots of flow cells at once. For example, ordering 48 MinION flow cells means that the per-flowcell price is \$500 USD instead of \$900 USD. That's 48 flow cells spread over the course of an entire year, with costs paid on delivery. If you can encourage multiple research groups to put in a combined bulk order, it's not too far out of the realm of possibility.
The long-read advantage
I've been told that a researcher is able to get 50M paired-end Illumina reads from a sequencing service centre for about $100 CAD.
50M paired end reads is about 10 Gb, which can be achieved on a MinION run (\$600 - \$1000), or about 1/10th of a PromethION run (\$850 - \$1500).
Even given those prices, it doesn't take much of a long-read advantage for MinION to become the most economical option. My own studies suggest that 1M full-length reads from nanopore sequencing (i.e. about 1/10th of a MinION R9.4.1 flow cell) gives results that are at least as good as 20M short paired-end Illumina reads.


text:
Finding the nearest gene to a specified gene region
I've got some data on genomic regions of interest in a bedgraph-like format from the Mus musculus [mm10] genome, and would like to find out the nearest gene to these regions.
Some of the regions will reside within a gene (in which case I'd like to know that they're within a gene), some regions may reside within multiple genes (in which case I'd like to know all gene names), and some regions are not within a gene (in which case I'd like to know the closest gene: the shortest sequence gap between either end of the gene and either end of the tagged genomic region. In this case I want to do an unstranded search, because it's possible that I might have regions describing suppressing transcripts. Here's some example data:
Ref         Start             End      fwd     rev
chrX     73716152        73716152        0      -2
chrX    167207094       167209162        0      -3
chrY     30770844        30772724        3       0
chr1     24613189        24613641        0       3
chr1     24613971        24614812        0       4
chr1     35788729        35789387       -5       0
chr1     93151586        93160225        0      -2
chr1    107511454       107525597       -2       0
chr1    130882067       130887413        0      -3
chr1    134182466       134189942        2       0
chr1    135818954       135832956       -2       0
chr1    149520498       149521157        0      -3
chr1    166166500       166166500       -2       0
chr2      5379059         5379374        0      -2
chr2     32379110        32381837        0      -2
chr2     32381854        32381858        0      -2
chr2     32384633        32384633        0      -2
chr2     32384636        32387773        0      -3
chr2    120539390       120563799        0      -2
chr2    125946785       125975527        2       0
chr2    127633311       127633684        0      -2
chr2    127634841       127656210        0      -3
chr2    164562717       164562717       -2       0
chr2    164562719       164562720       -2       0
chr2    164568498       164568499       -2       0
chr2    181188007       181188388       -2       0

I'd like to add another 'Gene' column to this table that gives the official names of the nearest gene(s) for each region defined by <Ref>:<Start>-<End>.

What would be the desired output format? Would you like the name of the gene in that new column? Are you interested in the nearest gene in the same strand or it doesn't matter? What happen if multiple genes are in those regions?

You're looking for bedtools closest, which thankfully defaults both considering overlapping genes and outputting all hits in case of ties.
  bedtools closest -a foo.txt -b genes.bed | awk '{foo = sprintf("%s:%i-%i", $6, $7, $8); OFS="\t"; print $1, $2, $3, $4, $5, foo}'

You'll need to strip the header and ensure the file is sorted. I've piped the output to awk to reformat things. The output then is something like:
chr1    24613189    24613641    0   3   chr1:24613188-24613971
chr1    24613971    24614812    0   4   chr1:24613973-24614651
chr1    35788729    35789387    -5  0   chr1:35806973-35810462
chr1    93151586    93160225    0   -2  chr1:93151348-93160948
chr1    93151586    93160225    0   -2  chr1:93151582-93160870
chr1    107511454   107525597   -2  0   chr1:107511422-107525598

Those results are for the Gencode M15 release (the BED file I used has transcripts rather than genes, so you'll get a bit different coordinates). As an aside, the gene name would normally be output as well, just add column 9 to the output.

With the prompting from @Devon_Ryan's answer, I came up with a solution that worked well for me. I downloaded the mm10 RefSeq gene locations from UCSC [Fields: chrom, txStart, txEnd, name2], and sorted both input files using bedtools sort. From there the bedtools closest command got me almost to where I wanted to go:
(echo -e "Chr\tStart\tEnd\tfwd\trev\tm_Chr\tm_Start\tm_End\tGene\tDistance";
  bedtools closest -d -a regions.bed -b NCBI_RefSeq_mm10.bed) > annotated.tsv

Output example:
Chr     Start   End     fwd     rev     m_Chr   m_Start m_End   Gene    Distance
chrX    73716152        73716152        0       -2      chrX    73686177        73716204        Bcap31  0
chrX    73716152        73716152        0       -2      chrX    73686177        73717815        Bcap31  0
chrX    73716152        73716152        0       -2      chrX    73686177        73716401        Bcap31  0
chrX    73716152        73716152        0       -2      chrX    73686177        73716502        Bcap31  0
chrX    167207094       167209162       0       -3      chrX    167207092       167213479       Tmsb4x  0
chrX    167207094       167209162       0       -3      chrX    167207093       167209218       Tmsb4x  0
chrY    30770844        30772724        3       0       chrY    30725674        30738546        Gm29582 32299
chr1    24613189        24613641        0       3       chr1    24257682        24587535        Col19a1 25655
chr1    24613189        24613641        0       3       chr1    24257682        24587535        Col19a1 25655
chr1    24613189        24613641        0       3       chr1    24257682        24587535        Col19a1 25655
chr1    24613189        24613641        0       3       chr1    24257682        24587535        Col19a1 25655
chr1    24613971        24614812        0       4       chr1    24257682        24587535        Col19a1 26437
chr1    24613971        24614812        0       4       chr1    24257682        24587535        Col19a1 26437
chr1    24613971        24614812        0       4       chr1    24257682        24587535        Col19a1 26437
chr1    24613971        24614812        0       4       chr1    24257682        24587535        Col19a1 26437
chr1    35788729        35789387        -5      0       chr1    35869594        35881164        1110002O04Rik   80208
chr1    93151586        93160225        0       -2      chr1    93151353        93160935        2310007B03Rik   0
chr1    93151586        93160225        0       -2      chr1    93151353        93160957        2310007B03Rik   0
chr1    93151586        93160225        0       -2      chr1    93151354        93160948        2310007B03Rik   0
chr1    93151586        93160225        0       -2      chr1    93151354        93160870        2310007B03Rik   0
chr1    107511454       107525597       -2      0       chr1    107500891       107525600       Serpinb2        0
chr1    107511454       107525597       -2      0       chr1    107511422       107525600       Serpinb2        0
chr1    107511454       107525597       -2      0       chr1    107511509       107525600       Serpinb2        0
chr1    130882067       130887413       0       -3      chr1    130882073       130887411       Il24    0
chr1    134182466       134189942       2       0       chr1    134182169       134190031       Chil1   0
chr1    134182466       134189942       2       0       chr1    134182403       134190031       Chil1   0
chr1    134182466       134189942       2       0       chr1    134182408       134190031       Chil1   0
chr1    135818954       135832956       -2      0       chr1    135818597       135833341       Lad1    0
chr1    149520498       149521157       0       -3      chr1    149410728       149449931       Gm29398 70568
chr1    166166500       166166500       -2      0       chr1    166130459       166166510       Gpa33   0
chr2    5379059         5379374         0       -2      chr2    5293456         5676046         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293456         5714762         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293456         5714762         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293456         5714762         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293458         5383533         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293458         5569754         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293458         5601302         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293458         5715350         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5298992         5715357         Camk1d  0
chr2    32379110        32381837        0       -2      chr2    32379100        32381915        1110008P14Rik   0
chr2    32381854        32381858        0       -2      chr2    32379100        32381915        1110008P14Rik   0
chr2    32384633        32384633        0       -2      chr2    32384636        32387739        Lcn2    3
chr2    32384636        32387773        0       -3      chr2    32384636        32387739        Lcn2    0
chr2    120539390       120563799       0       -2      chr2    120506829       120563831       Zfp106  0
chr2    120539390       120563799       0       -2      chr2    120507433       120563851       Zfp106  0
chr2    120539390       120563799       0       -2      chr2    120507433       120563853       Zfp106  0
chr2    125946785       125975527       2       0       chr2    125859108       125984298       Galk2   0
chr2    125946785       125975527       2       0       chr2    125859228       125983587       Galk2   0
chr2    125946785       125975527       2       0       chr2    125866111       125983587       Galk2   0
chr2    125946785       125975527       2       0       chr2    125866218       125983587       Galk2   0
chr2    125946785       125975527       2       0       chr2    125866219       125983587       Galk2   0
chr2    125946785       125975527       2       0       chr2    125866222       125984298       Galk2   0
chr2    125946785       125975527       2       0       chr2    125920565       125983587       Galk2   0
chr2    127633311       127633684       0       -2      chr2    127633225       127656695       Mal     0
chr2    127633311       127633684       0       -2      chr2    127633225       127656695       Mal     0
chr2    127634841       127656210       0       -3      chr2    127633225       127656695       Mal     0
chr2    127634841       127656210       0       -3      chr2    127633225       127656695       Mal     0
chr2    164562717       164562717       -2      0       chr2    164562560       164568510       Wfdc2   0
chr2    164562717       164562717       -2      0       chr2    164562715       164568506       Wfdc2   0
chr2    164562719       164562720       -2      0       chr2    164562560       164568510       Wfdc2   0
chr2    164562719       164562720       -2      0       chr2    164562715       164568506       Wfdc2   0
chr2    164568498       164568499       -2      0       chr2    164562560       164568510       Wfdc2   0
chr2    164568498       164568499       -2      0       chr2    164562715       164568506       Wfdc2   0
chr2    181188007       181188388       -2      0       chr2    181187221       181188506       Ppdpf   0
chr2    181188007       181188388       -2      0       chr2    181187342       181188504       Ppdpf   0

After that it was just a matter of aggregating the gene names for each unique region with a quick R script:
library(dplyr);
data.tbl <-
    read.delim("annotated.tsv") %>%
    group_by(Chr, Start, End, fwd, rev) %>%
    summarise(Genes = paste(unique(Gene), collapse=","), Distance=min(Distance));

write.table(data.tbl, "collapsed_annotated.tsv",
            sep="\t", row.names=FALSE, quote=FALSE);

Collapsed output example:
Chr     Start   End     fwd     rev     Genes   Distance
chr1    24613189        24613641        0       3       Col19a1 25655
chr1    24613971        24614812        0       4       Col19a1 26437
chr1    35788729        35789387        -5      0       1110002O04Rik   80208
chr1    93151586        93160225        0       -2      2310007B03Rik   0
chr1    107511454       107525597       -2      0       Serpinb2        0
chr1    130882067       130887413       0       -3      Il24    0
chr1    134182466       134189942       2       0       Chil1   0
chr1    135818954       135832956       -2      0       Lad1    0
chr1    149520498       149521157       0       -3      Gm29398 70568
chr1    166166500       166166500       -2      0       Gpa33   0
chr2    5379059         5379374         0       -2      Camk1d  0
chr2    32379110        32381837        0       -2      1110008P14Rik   0
chr2    32381854        32381858        0       -2      1110008P14Rik   0
chr2    32384633        32384633        0       -2      Lcn2    3
chr2    32384636        32387773        0       -3      Lcn2    0
chr2    120539390       120563799       0       -2      Zfp106  0
chr2    125946785       125975527       2       0       Galk2   0
chr2    127633311       127633684       0       -2      Mal     0
chr2    127634841       127656210       0       -3      Mal     0
chr2    164562717       164562717       -2      0       Wfdc2   0
chr2    164562719       164562720       -2      0       Wfdc2   0
chr2    164568498       164568499       -2      0       Wfdc2   0
chr2    181188007       181188388       -2      0       Ppdpf   0
chrX    73716152        73716152        0       -2      Bcap31  0
chrX    167207094       167209162       0       -3      Tmsb4x  0
chrY    30770844        30772724        3       0       Gm29582 32299

Another option is BEDOPS closest-features, which efficiently locates nearest annotations, overlapping or otherwise.
$ closest-features --closest --delim '\t' intervals.bed genes.bed > answer.bed

It's easy to reformat output by piping to awk, say.
If files are very large, you can parallelize work into "map-reduce" form by specifying the chromosome name as an option:
$ closest-features --chrom chr1 --closest --delim '\t' intervals.bed genes.bed > answer.chr1.bed
$ closest-features --chrom chr2 --closest --delim '\t' intervals.bed genes.bed > answer.chr2.bed
...

To reduce, just bedops -u all the answer.*.bed files.
Running these per-chromosome jobs in parallel makes the work take about as long as searching chr1, which tends to be the largest chromosome.


text:
Removing repeated reads from nanopore 1D² reads
This is summary of discussions and a question that was posted by Devon O'Rourke on Twitter
Following Albacore basecalls on a 1D² library I get two sets of .fq files, summary stats, etc.: one for the 1D basecalling script, and one for the 1D² script. I'd like to use the reads generated from these scripts in a genome assembly.
After a little bit of grep searching, it seems like there are overlapping reads among the 1D and 1D² directories - namely, those 1D reads that generated 1D² reads!
How do people approach the choice of selecting which directory to pull the unique read from if you're planning on assembling a genome? I'd think that the 1D² reads are of higher quality and be preferred, but perhaps there are 1D² reads which are lower in quality than one of the two 1D pairs?
Just looking for advice on how folks approach 1D² libraries, fastq outputs, and subsequent mapping. For what it's worth I'm hoping to use Minimap2 to proceed with alignments but am happy to hear others thoughts on alternatives.
Update
I tried the first solution proposed by @gringer, but my filtered_1d_reads.fq file contains same number of reads as raw_1d_reads.fq file.
I'm pretty sure the problem is the naming convention difference between the linear and paired reads. linear reads always follow a read_id of alphanumeric (8)-(4)-(4)-(4)-(12), while paired reads are (8)-(4)-(4)-(4)-(20)-(4)-(4)-(4)-(12). When I looked at it in more detail, the paired read "read_id" naming is just a concatenation of the two linear "read_id"'s that are partnered. So that (20) value is really the last (12) of the first linear read, and the first (8) of the second linear pair.

Bit silly, but here it goes. > I'd think that the 1D² reads are of higher quality and be preferred,

but perhaps there are 1D² reads which are lower in quality than one of
the two 1D pairs? Based on this it would make the most sense to have a script preferentially select 1D2 reads unless they have a lower quality than the 1D 'parents', and keep all 1D reads that did not contribute to a 1D2 read. I don't have a dataset to play with, but this shouldn't be too hard to write I believe.

An additional QC step (as suggested by @Wouter_De_Coster) would be a good idea. It could probably be done by using the run statistics file in the workspace directory to generate the read ids for filtering.

Parsing the sequencing_summary.txt is convenient and fast. I have an old 2D example here in my test data folder, showing that it contains qscores of template, complement, and 2d. If that would be the same for 1D then filtering out 1D2 reads which we don't want to keep because they have lower quality than parents would be a breeze.

You can exclude the 1D² reads by generate a list of the ids from the 1D² reads, then filtering those reads out from the 1D fastq files. Here's one way to do that using some accessory scripts that I have written:
$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | awk '{print $2}' > readNames.txt
$ cat 1d_reads.fq | ~/scripts/fastx-fetch.pl -v -i readNames.txt > filtered_1d_reads.fq

In order to deal with the concatenated read IDs, it's necessary to switch to a language like Perl that supports the use of backreferences for the splitting. The following perl command (replacing the first line in the above procedure) inserts a line break after the 36th character in the 2nd value, then prints that value:
$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | perl -lane '$F[1] =~ s/(.{36})/$1\n/; print $F[1];' > readNames.txt

So a full sequence of operations would look something like this:
$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | perl -lane '$F[1] =~ s/(.{36})/$1\n/; print $F[1];' > readNames.txt
$ cat 1d_reads.fq | ~/scripts/fastx-fetch.pl -v -i readNames.txt > filtered_1d_reads.fq
$ cat 1d2_reads.fq filtered_1d_reads.fq > unique_reads.fq


text:
How can I improve a long-read assembly with a repetitive genome?
I'm currently trying to assembly a genome from a rodent parasite, Nippostrongylus brasiliensis. This genome does have an existing reference genome, but it is highly fragmented. Here are some continuity statistics for the scaffolds of the current Nippo reference genome (assembled from Illumina reads):
Total sequences: 29375
Total length: 294.400206 Mb
Longest sequence: 394.171 kb
Shortest sequence: 500 b
Mean Length: 10.022 kb
Median Length: 2.682 kb
N50: 2024 sequences; L50: 33.527 kb
N90: 11638 sequences; L90: 4.263 kb

This genome is most likely difficult to assemble because of the highly repetitive nature of the genomic sequences. These repetitive sequences come in (at least) three classes:

Tandem repeats with a repeat unit length greater than the read length of Illumina sequencers (e.g. 171bp)
Tandem repeats with a cumulative length greater than the fragment length of Illumina sequencers, or the template length for linked reads (e.g. 20kb)
Complex (i.e. non-repetitive) sequence that appears at multiple places throughout the genome

Canu seems to deal quite well with the first two types of repeats, despite the abundance of repetitive structure in the genome. Here's the unitigging summary produced by Canu on one of the assemblies I've attempted. Notice that about 30% of the reads either span or contain a long repeat:
category            reads     %          read length        feature size or coverage  analysis
----------------  -------  -------  ----------------------  ------------------------  --------------------
middle-missing        694    0.07     7470.92 +- 5552.00        953.06 +- 1339.13    (bad trimming)
middle-hump           549    0.05     3770.05 +- 3346.10         74.23 +- 209.86     (bad trimming)
no-5-prime           3422    0.33     6711.32 +- 5411.26         70.92 +- 272.99     (bad trimming)
no-3-prime           3161    0.30     6701.35 +- 5739.86         87.41 +- 329.42     (bad trimming)

low-coverage        27158    2.59     3222.51 +- 1936.79          4.99 +- 1.79       (easy to assemble, potential for lower quality consensus)
unique             636875   60.76     6240.20 +- 3908.44         25.22 +- 8.49       (easy to assemble, perfect, yay)
repeat-cont         48398    4.62     4099.55 +- 3002.72        335.54 +- 451.43     (potential for consensus errors, no impact on assembly)
repeat-dove           135    0.01    16996.33 +- 6860.08        397.37 +- 319.52     (hard to assemble, likely won't assemble correctly or even at all)

span-repeat        137927   13.16     9329.94 +- 6906.27       2630.06 +- 3539.53    (read spans a large repeat, usually easy to assemble)
uniq-repeat-cont   155725   14.86     6529.83 +- 3463.16                             (should be uniquely placed, low potential for consensus errors, no impact on assembly)
uniq-repeat-dove    28248    2.70    12499.99 +- 8446.95                             (will end contigs, potential to misassemble)
uniq-anchor          5721    0.55     8379.86 +- 4575.71       3166.22 +- 3858.35    (repeat read, with unique section, probable bad read)

However, the third type of repeat is giving me a bit of grief. Using the above assembly, here are the continuity parameters from the assembled contigs:
Total sequences: 3505
Total length: 322.867456 Mb
Longest sequence: 1.762243 Mb
Shortest sequence: 2.606 kb
Mean Length: 92.116 kb
Median Length: 42.667 kb
N50: 417 sequences; L50: 194.126 kb
N90: 1996 sequences; L90: 35.634 kb

It's not a bad assembly, particularly given the complexity of the genome, but I feel like it could be improved by tackling the complex genomic repeats in some fashion. About 60Mb of the contigs in this assembly are linked with each other in a huge web (based on the GFA output from Canu):

The repetitive regions are typically over 500bp in length, average about 3kb, and I've seen at least one case which seems to be a 20kb sequence duplicated in multiple regions.
The Canu defaults seem to give the best assembly results for the few parameters that I've tried, with one exception: trimming. I've tried playing around a little bit with the trimming parameters, and curiously a trimming coverage of 5X (with overlap of 500bp) seems to give a more contiguous assembly than with a trimming coverage of 2X (with the same overlap).
If anyone is interested in having a look at these data themselves, called FASTQ files from Nippo sequencing runs can be found here. Raw nanopore signal files are available within ENA project PRJEB20824. There's also a Zenodo archive here that contains the GFA and assembly contigs.
I can use Illumina data to correct the Canu assembly, but that doesn't help with resolving the "type 3" repeats. The regions are sufficiently similar that illumina reads get mapped to multiple points in the genome. The Illumina contigs are high quality (i.e. they have good BUSCO scores, indicating few variant errors), but quite short. Any sniff of a repeat and the contig ends. I've got more than a few examples of regions that would make an Illumina read (even 10x linked reads) cower in fear.
Does anyone have any other suggestions on how I could resolve these complex repeats?
Computational solutions would be preferred, but resequencing is not out of the question.

"A few" 100kb reads won't help much. You need to apply the ultra-long protocol, which is different from the standard protocol.
You can't resolve 20kb near identical repeats/segdups with 10kb reads. All you can do is to bet your luck on a few excessively long reads spanning some units by chance. For divergent copies, it is worth looking at this paper. It uses Illumina reads to identify k-mers in unique regions and ignores non-unique k-mers at the overlapping stage. The paper said that this strategy is better than using standard overlappers, which I buy, but probably it can't resolve a 20kb segdup with a handful of mismatches, either.
Such mismatch-based approaches always have limitations and may not work for recent segdups/repeats. The ultimate solution is to get long reads, longer than your repeat/segdup units. The ~100kb reads in the recent preprint will be a game changer for you. If your ~20kb repeats are not tandem, 10X's ~100kb linked reads may help, too.


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

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

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

For us, our best cDNA sequencing run on MinION produced over 15M reads in June 2020. We've since shifted over to doing cDNA sequencing on PromethION flow cell (using a P2 Solo), and our best so far with that run was 89M reads.

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

The DNA cleanup
The library preps

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

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


text:
What is the difference between NGS, 2GS, SBS and HTS?
I've come across a bit of confusion about the initialism NGS, so think it would be a good idea to clarify this term (and similar terms like 2GS, SBS, and HTS) for this site with a bit of discussion.
What are the most common initialisms used to describe different sequencing technologies?

I think your answer is good, but also feel these terms are mainly used for marketing; much more informative to name the sequencing company and platform. NGS is the worst, some 'NGS' like SOLID and 454 are basically retired

NGS isn't as bad as SMRT, which, despite the expansion of "single-molecule real-time sequencing", is considered specific to PacBio and frowned upon as a tag for nanopore sequencing (which is also single-molecule and real-time).

I am not sure how appropriate it is at this point to still refer to sequencing as next-generation sequencing. The leading NGS technology is Illumina/Solexa that has been around for over 10 years at this point. 454 was around even earlier. It's not really "next" at this point.
Opinions aside, I would also add "third-generation sequencing" to that list, referring to long read technologies like PacBio and Oxford Nanopore. See this question on Biology SE for more details.

I have a bit of trouble classifying PacBio in a sequencing generation. It is a model-based sequencing technology (dependent on a polymerase and pre-mixed nucleotides for sequencing), but produces single-molecule reads sequentially over the course of a run. I therefore consider it to be somewhere between second and third generation sequencing technologies.

Here are my attempts at definitions:
Sanger: A method of sequencing that depends on chain-terminatiing dideoxynucleotides. This sequencing uses the differential flow of DNA sequences of different lengths through a gel to determine the original DNA sequence, producing a single sequence per reaction container.
NGS: Next-generation sequencing, also referred to as 2GS (second-generation sequencing). This term is used to describe the first wave of sequencing technologies that followed Sanger sequencing technology. The use of NGS has become more confusing with the advent of long-read sequencing, because it's a common assumption that "next-generation" refers to the most recent technology (which is incorrect in this case).
HTS: High-throughput sequencing. This term describes any type of sequencing technology that produces large amounts of data, usually in the form of millions of different sequences produced from the same sequencing run.
SBS: Sequencing by synthesis. This term describes a method of sequencing that depends on the synthesis of [DNA] bases in order for sequencing to be carried out. This definition can extend into long-read sequencing (e.g. PacBio sequencers depend on synthesis during sequencing), but is more typically associated with only the second-generation sequencing technology.

For me ngs and hts are the same thing - any massively parallel sequencing (anything after sanger). That is further divided to 2nd generations sequencing (short reads) and 3rd generation sequencing (long reads). Sequencing by synthesis can by whatever sequencing that uses template and synthesis of free nucleotides (454, Illumina, PacBio).

If you have to divide "next-generation" into multiple generations, it's probably not a good definition.

Some people will distinguish the ABI SOLiD and Nanostring Hyb+Seq sequencing technologies as "sequencing by ligation" rather than "sequencing by synthesis". I currently lump these all together as SBS because I haven't found a better descriptive name for the collective term.

I did bit of digging about "Sequencing by synthesis" : it is defined on an Illumina flyer, page 15, it is defined though florescent died nucleotides... PacBio write on their web : "Also, while it is a form of ‘sequencing by synthesis’, it operates in real-time rather than incorporating a single base (or type of base) at a time."

It's not quite "real-time" sequencing, as the base-calling is done by decoding fluorescence "movies" for each photodetector (ZMWG)... but then again nothing is truly real-time unless it's detecting and calling single bases at the same time.

Also, look at answer on very similar question on biology.SE. What you write here is really not what the most of scientist see behind these terms...

The terms "second-generation sequencing", "third-generation sequencing", "short-read sequencing", and "long-read sequencing" are used commonly but (as far as I know) lack commonly used acronyms/initialisms.
In the forensic genetics community, the acronym MPS is used commonly to refer to "massively parallel sequencing." This definitely includes 2G/short-read sequencing (such as Illumina). It's less clear whether MPS includes 3G/long-read sequencing such as PacBio or Oxford Nanopore, but since these are rarely used in a forensics context it's a moot point.
Of course, it's important to acknowledge (as @burger has done) that there is a lot of imprecision and opinion regarding the appropriate scope of each term. Some folks seem to want a term to refer to everything post-Sanger, which includes an incredible diversity of platforms with respect to sequencing technology, read length, throughput, and error profiles. Some very reasonable questions should make us think seriously about our chosen terminology: How long do (for example) Illumina reads have to get before they're no longer classified as short reads? How much throughput is required of a sequencer to be classified as "high" throughput? When do these distinctions really matter?

authors: ลูกชิ้นเด้ง เตี๋ยวหนับน้ำตกตับ


text:
How do I increase the sensitivity of Salmon Alevin?
I'm in a little bit of a bind with targeted single-cell sequencing. I'm trying to match up our reads to the targeted amplicon panel (418 targets), and all but one have matched successfully with Salmon Alevin. I have a transcript that produces many 75bp reads that only match the first 28bp of a target, and Alevin is ignoring those reads as a match (the 28bp sequence is genome-unique).
Here's an example read, with matched sequence:
qual FFFFFFFFFFFFFFFFFFFFFFFFFFFF,,,::,,:,:,,,,,,,,,,,,,,,,,:,,,,,,,,:,,,,,:,:
read GCTGCATAATCTCTTCATTCCGAGGAGCACCCCACGCGCGCCGCGCCCCCCCGCCCCCGGGGGGGGGGGGGGG
SNPs ...................................G.G.G..G.GC.....CG...CCGGGGGGGGG..GG..
 ref GCTGCATAATCTCTTCATTCCGAGGAGCACCCCACCCCCACCCCCACCCCCACCCCATTCTTAAATTGTTTGG

As you may be able to see from the above sequence, the reads for this particular transcript get up to 28bp (sometimes with one or two errors), then the wheels start coming off the sequencer, showing an abundance of Gs. This was a NovaSeq run, and a polyG read is the same as no signal using its two-dye process, so I suspect that the sequencing stops after ~35bp, rather than actually synthesising Gs from thin air.
I have lots of these reads; there are about 400k reads that match the reference transcript (Cxcr4, just in case it matters), and about 250k of those reads have a G homopolymer that's at least 8bp in length.
Any ideas on what to do?
FWIW, the above match was taken from a Bowtie2 alignment I did on the reads. I could use Alevin for cell barcode correction only and do the UMI deduplication manually, but it seems a shame to do that when it's only one transcript (from 418) that's spoiling a fully-automated process.
I did try adding in an additional transcript that included the full sequence of a short-match read, but there seem to be too many errors in the non-matching sequence for any reads to count (at least... I assume that's why the counts are not coming through; I'm not quite sure how I can work out what's going on behind the scenes).
Is there any way to tweak Salmon Alevin to get it to accept these shorter matches? If I can change the kmer size, how do I do that?
End trimming is unlikely to help, because that's going to shorten the reads and make them less likely to map to the target sequence. I want a way to include these reads as valid hits to the target sequence.

So, I think there are a few potential options here.  Alevin is using selective-alignment internally to determine the mappings for the reads.  So, even if you have a k-mer supporting the mapping, if you have gibberish for the rest of the read, the mapping score is going to be very poor and you're not going to recover that mapping locus.  There are some parameters you can mess with to try to increase the sensitivity of mapping, though you should always be aware that upping the sensitivity a lot will mean there is less aggressive filtering early in the mapping process which can slow things down.  The main parameter you should tweak to allow lower alignment scores is --minScoreFraction; the smaller this value the lower (worse) will be the permitted scores for valid mappings.  Additionally, you can combine changes to the --minScoreFraction with the --softclip option to simply score the non-matching end of read as 0 (rather than a long deletion or large set of mismatches); this will lower the effect of bad read ends.
Another alternative you might consider is to try out alevin-fry for processing this data, and mapping the reads with the --sketch flag.  This will only check certain structural constraints of the mapping and won't force a particular minimum alignment score.  The caveat here is that dropping the requirement for alignment scoring can increase the rate of spurious mappings (reads deriving from elsewhere that spuriously map to some indexed target simply by virtue of sharing some sequence).  In transcriptome wide analyses, this can be greatly mitigated by building an appropriate index and using a UMI resolution strategy that complements this expanded index. We've not examined this in detail for targeted single-cell sequencing.  However, if you want to give it a try and run into any issues, feel free to drop an issue over on GitHub and we'd be happy to help work through it.

Thank you; tweaking the minScoreFraction did the trick. I went from a max Cxcr4 read count per cell of 1 with default parameters, to 4 with minScoreFraction=0.5, to 26 with minScoreFraction=0.33.

@ATPoint commented that the default kmer size in Salmon indexing is 31, which led me to a solution that involved altering the index generation step:
Instead of:
salmon index -t merged_salmon_${targetName}.fasta -i salmonIndex_${targetName}

I now have:
salmon index -t merged_salmon_${targetName}.fasta -i salmonIndex_k11_${targetName} -k 11

In other words, I changed the kmer size from its default of 31 to 11, which matches the BD-recommended seed substring length for the standard Bowtie2 approach ("-L 10", but I increased that to 11 because Salmon needs odd kmer sizes). I'm pretty sure that my 28bp-mapping reads won't be picked up with a kmer size of 31. Because this is targeted sequencing, I don't expect that a reduced kmer size will be an issue, but I'll keep a watch out for that.
I'm currently re-running the analysis. If it still doesn't work, then I'll see about adding in the tweaks suggested by @nomad.
Update: still not working, so I'll shift onto score tweaks.


text:
Finding the location and unit length of repetitive sequences within a long read
After discovering a few difficulties with genome assembly, I've taken an interest in finding and categorising repetitive DNA sequences, such as this one from Nippostrongylus brasiliensis [each base is colour-coded as A: green; C: blue; G: yellow; T: red]:

[FASTA file associated with this sequence can be found here]
These sequences with large repeat unit sizes are only detectable (and assembleable) using long reads (e.g. PacBio, nanopore) because any subsequence smaller than the unit length will not be able to distinguish between sequencing error and hitting a different location within the repeat structure. I have been tracking these sequences down in a bulk fashion by two methods:

Running an all-vs-all mapping, and looking for sequences that map to themselves lots of times
Carrying out a compression of the sequence (e.g. bzip2), and finding sequences that have a compression rate that is substantially higher than normal

After I've found a suspicious sequence, I then want to be able to categorise the repeat (e.g. major repeat length, number of tandem repeats, repetitive sequence). This is where I'm getting stuck.
For doing a "look, shiny" demonstration, I currently have a very manual process of getting these sequences into a format that I can visualise. My process is as follows:

Use LAST to produce a dot plot of self-mapping for the mapping
Visually identify the repetitive region, and extract out the region from the sequence
Use a combination of fold -w <width> and less -S to visually inspect the sequence with various potential repeat unit widths to find the most likely repeat unit size
Display the sequence in a rectangular and circular fashion using my own script, wrapping at the repeat unit length

But that process is by no means feasible when I've got thousands of potential repetitive sequences to fish through.
Is there any better way to do this? Given an arbitrary DNA sequence of length >10kb, how can I (in an automated fashion) find both the location of the repeat region, and also the unit length (bearing in mind that there might be multiple repeat structures, with unit lengths from 30bp to 10kb)?
An example sequence can be found here, which has a ~21kb repeat region with ~171bp repeat units about 1/3 of the way into the sequence.
A Kmer-based Analysis
I've now seen human sequences with repetitive regions in excess of 10kb (i.e. out of the range of most linked-read applications). My current idea is centred around creating hash tables of short sequences (currently 13-mers) and tracking their location:

Process the sequence, storing the location of each kmer
For each kmer, find out how many times it appears in the sequence
For repeated kmers, find out how much of a gap there is between the next time that kmer occurs
Report the median and modal gap length of repeated kmers, and statistics associated with their frequency in the sequence

Some local repetitive regions may be lost in the statistics with this approach, it's hard to tell if there are multiple repetitive regions within a single sequence, and if the repeat units are themselves slightly repetitive (enough that a kmer is duplicated within a repeat unit), then the algorithm will under-report the repetitiveness (see step 3).

Would RepeatMasker help you here? Or does it only find known repetitive elements and regions of low complexity? I've never tried to use it to identify simple repeats.

These repeats have high-complexity repeat units; I haven't found the regions being marked by tantan or dust, and suspect that RepeatMasker will be similar. The web service at least isn't able to find the repeat region in the test sequence that I've linked.

I know about this blog post from Gene Myers where he is using strategy of masking repeats in long reads to speed up genome assembly. I am not sure if it is exactly what you want, but I thought it might be useful.

Thank you. TANmask (at the very end of the post) does look very much like the sort of algorithm I'm looking for. It will at least solve the boundary problem, if not the repeat unit length statistic. Can you write that up as an answer?

I've sorted the visualisation out. Here are three alternative representations of repetitive structures for the same sequence:

These were generated using the same R script, callable from the command line:
 $ repaver.r -style dotplot -prefix dotplot MAOA_chrX_43645715-43756016.fa
 $ repaver.r -style profile -prefix profile MAOA_chrX_43645715-43756016.fa
 $ repaver.r -style semicircular -prefix semi MAOA_chrX_43645715-43756016.fa

More details about this are in the presentation I gave at Queenstown Research Week, 2018. I also wrote a chapter in a peer-reviewed eBook here.
This is fast enough that I can run it on the nanopore C. elegans genome in about half an hour, producing these plots for each contig. I don't quite have a method to iterate through this plot and pick out the dominant repeat length at each location, but that's a [relatively] simple extension of what I've already done.
With a lot of optimisations for speed and memory consumption, I've now been able to run it on the full human genome. It takes a couple of days on my desktop computer (64GB RAM + SSD swap) to categorise the 100bp repeat structure of the assembled T2T/CHM13-v1.0 chromosomes.
Code here.

It could be an idea to fragment the long reads into small sequences, like simulating Illumina reads of 150 bp, and then map these small sequences against the original long reads and extract regions with a high coverage?


text:
Can second search with limited search space improve result in MS proteomics data?
First, please know that I am a novice in this field and that it is thus very likely that you will regard this question as basic or flat stupid. Please be gentle and don't be afraid to criticize. I am asking this so I can learn. Also, don't feel afraid to discuss in comments if you can't give a straight answer; all input is useful.
I am currently working in a project where we have fractionated samples in order to dig deeper into the proteome. The fractionated samples have been analysed with shotgun LC-MS and then searched with a contemporary software using a whole proteome database and to the reversed data base in order to estimate FDR. The search result is cut-off so that the highest scoring sequences are included with an estimated FDR of 1%. In this project we have 150 samples in which we can expect to find certain proteins in almost all samples, while others are unlikely to be found in any of them (but we don't know which yet). If we were to decrease the search space to include only relevant proteins (which as previously stated are not known to us yet), we would boost our rate of true positives included after FDR cut-off, as we decrease the search space which noise randomly might match to with a high score (right?). As we can't know which ones are relevant beforehand, would it be possible to determine relevant proteins with a preliminary search? That is, to do as follows:

Search against whole proteome with a certain FDR cut-off (possibly larger than 1%)
Make a new data base including only proteins that has been found in any of the samples.
Search data against new data base.

To me it seems like there is a possibility that this could increase the number of true positives that are found, but I also feel like I might be missing something crucial. Please, enlighten me: Is this methododology unbiased and valid? What risks may there be? Anything more one has to think about before attempting something like this?
The LC-MS was a top 10 shotgun MSMS using a quadrupole-Orbitrap instrument. The search was done using MaxQuant where fragmentation spectra are use to identify MS1 peaks, which intensity is registered as a quantitative datum. All files were searched together but are configured to be separated in the results file in regards to both sample and fraction. We will likely not perform additional analysis of these samples.


text:
How can longest isoforms (per gene) be extracted from a FASTA file?
Is there a convenient way to extract the longest isoforms from a transcriptome fasta file? I had found some scripts on biostars but none are functional and I'm having difficulty getting them to work.
I'm aware that the longest isoforms aren't necessarily 'the best' but it will suit my purposes.
The fasta was generated via Augustus. Here is what the fasta file looks like currently (sequence shortened to save space)
>Doug_NoIndex_L005_R1_001_contig_2.g7.t1
atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgattcgtatgatcaactcatca
ttaatataaccaaaaacctagaaattctagccttcgatgatgttgcagctgcggttcttgaagaagaaagtcggcgcaagaacaaagaagatagaccg
>Doug_NoIndex_L005_R1_001_contig_2.g7.t2
atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgattcgtatgatcaactcatca

The format is as such:
Gene 1 isoform 1
Gene 1 isoform 2
Gene 2 isoform 1
Gene 2 isoform 2

and so forth. There are several genes that have more than one pair of isoforms (up to 3 or 4). There are roughly 80,000 total transcripts, probably 25,000 genes. I would like to extract the single longest isoform for each gene.

It depends on the exact format of your fasta. Where did you get it? What do the fasta header lines look like? Can you paste a small example?

Doug_NoIndex_L005_R1_001_contig_4.g13.t1

(sequence data here)

Doug_NoIndex_L005_R1_001_contig_4.g13.t2

(sequence data here)

What exactly have you tried? 2. What does them not being functional mean? 3. Do you just want the single longest one or the top N longest?

I've tried all of the solutions here: https://www.biostars.org/p/107759/

Non functional in that I couldn't get it work

I would like the single longest read for each set of isoforms

@ZincFingers: Please put all of this in your original post and include whether isoforms are ever separated from each other by unrelated genes.

You don't have the same sequence IDs there. Should we assume that the part after the last . in the ID line needs to be discarded?

I'm not sure what you mean. The ".t1" or ".t2" indicates which isoform the read is from, I would like to discard the shorterst isoforms and retain the longest isoform.

OK, so they do need to be removed.

While the solution from https://bioinformatics.stackexchange.com/users/96/daniel-standage should work (after adjusting for possible python3 incompatibility), the following is a shorter and less memory demanding method that uses biopython:
#!/usr/bin/env python
from Bio import SeqIO
import sys

lastGene = None
longest = (None, None)
for rec in SeqIO.parse(sys.argv[1], "fasta"):
    gene = ".".join(rec.id.split(".")[:-1])
    l = len(rec)
    if lastGene is not None:
        if gene == lastGene:
            if longest[0] < l:
                longest = (l, rec)
        else:
            lastGene = gene
            SeqIO.write(longest[1], sys.stdout, "fasta")
            longest = (l, rec)
    else:
        lastGene = gene
        longest = (l, rec)
SeqIO.write(longest[1], sys.stdout, "fasta")

If you saved this as filter.py, then filter.py original.fa > subset.fa would be the command to use.

Here is a solution in R. Could get really slow with big files. Works for the example you posted.
library(Biostrings)

## read your fasta in as Biostrings object
fasta.s <- readDNAStringSet("sample.fa")

## get the read names (in your case it has the isoform info)
names.fasta <- names(fasta.s)

## extract only the relevant gene and isoform id (split name by the period symbol)
gene.iso <- sapply(names.fasta,function(j) cbind(unlist(strsplit(j,'\\.'))[2:3]))

## convert to good data.frame = transpose result from previous step and add relevant column names
gene.iso.df <- data.frame(t(gene.iso))
colnames(gene.iso.df) <- c('gene','isoform')

## and length of isoforms
gene.iso.df$width <- width(fasta.s)

## split data.frame into list with entry for each gene
gene.iso.df.split <- split(gene.iso.df,gene.iso.df$gene)

## optional to keep all the information but really you just need indices
##gene.iso.df.split.best <- lapply(gene.iso.df.split,function(x) x[order(x$width)[1],])

## pull out the longest isoform ID for each gene (in case of a tie just take the first one)
best.id <- sapply(gene.iso.df.split,function(x) row.names(x)[order(x$width)[1]])

## subset your original reads with the subset
fasta.s.best <- fasta.s[best.id]

## export new fastafile containing longest isoform per gene
writeXStringSet(fasta.s, filepath='sample_best_isoform.fasta')

Late to the party here, but I like to try to avoid writing scripts when some command line magic will do. It's good practice to index your FASTA so use it.
samtools faidx <myfasta.fa>

A little bit of awk and sort can determine the largest isoform of each gene (and they don't even have to be sorted by name in the source FASTA).
awk -F'[\t.]' '{print $1,$2,$3,$4}' <myfasta.fa>.fai | sort -k4nr,4 | sort -uk1,2 | cut -f1-3 -d' '| tr ' ' '.' > selection.ls

The -F[\t.] specifies both the fai tab delimiters, and the full stops in your contig names (Doug_NoIndex_L005_R1_001_contig_2.g7.t1) as delimiters for awk.
For each line we print four columns. $1 to $3 are the three components of your contig name (after splitting on the .), and column $4 is the column in the index for sequence length.
We use sort to organise the lines by $4 (the contig length). Reverse sort with -r so the longest genes appear first.
We then subsort by columns $1 and $2 (the contig and gene name (g1...)). -u then drops any repeat of column $1 and $2 pairs (that is, all lines after the first pair -- the longest length).
Finally, we reassemble $1 to $3 with cut and tr to give you a list of sequence names.

A loop with our old friend samtools faidx can put your chosen sequences out to a new FASTA.
while read contig;
    do samtools faidx <myfasta.fa> $contig >> selection.fa;
done < selection.ls

Don't forget to index it, they are useful!
samtools faidx selection.fa

Oh and if you wanted to remove the .tN part, just switch out that cut -f1-3 -d' ' for cut -f1-2 -d' '

This solution is will work for your example and probably all Augustus-derived Fastas, but mileage will vary beyond that.
#!/usr/bin/env python
from __future__ import print_function
import sys

def parse_fasta(data):
    """Stolen shamelessly from http://stackoverflow.com/a/7655072/459780."""
    name, seq = None, []
    for line in data:
        line = line.rstrip()
        if line.startswith('>'):
            if name:
                yield (name, ''.join(seq))
            name, seq = line, []
        else:
            seq.append(line)
    if name:
        yield (name, ''.join(seq))

isoforms = dict()
for defline, sequence in parse_fasta(sys.stdin):
    geneid = '.'.join(defline[1:].split('.')[:-1])
    if geneid in isoforms:
        otherdefline, othersequence = isoforms[geneid]
        if len(sequence) > len(othersequence):
            isoforms[geneid] = (defline, sequence)
    else:
        isoforms[geneid] = (defline, sequence)

for defline, sequence in isoforms.values():
    print(defline, sequence, sep='\n')

Save as longest.py, and invoke on the command line like so.
python longest.py < intput.fasta > output.fasta

The line geneid = '.'.join(defline[1:].split('.')[:-1]) is key. Let me break it down.

defline[1:]: ignore the first character of the defline (the > symbol)
defline[1:].split('.'): split the string on . symbol
defline[1:].split('.')[:-1]: ignore the last value after the split (the isoform name)
'.'.join(defline[1:].split('.')[:-1]): join the split values again

This value will be the same for all isoforms for the same gene, so we use it in the script to keep track of the longest sequence associated with each gene ID.

This could save a lot of memory and a bit of time if you write out and reset isoforms whenever geneid not in isoforms. OP says isoforms are always next to each other.

Hi Daniel, I tried running this script and received the following error message:
File "longest.py", line 29
     print(defline, sequence, sep='\n')
                            ^
SyntaxError: invalid syntax

@ZincFingers: Change it to print(defline) and print(sequence) on separate lines.

Python 2 / 3 incompatibility. Will update.

An ex-coworker (Josep Avril) has written a couple of very useful little scripts that convert fasta to tbl (seqID<TAB>Sequence) and back again. These are extremely handy for this sort of thing (and are included at the end of this answer). Using them, you can convert your fasta to a one sequence per line format, keep the longest sequence with a simple awk script and convert back to fasta.
I am assuming that the part after the last . in the ID line should be removed (so that Doug_NoIndex_L005_R1_001_contig_2.g7.t1 and Doug_NoIndex_L005_R1_001_contig_2.g7.t2 both map to Doug_NoIndex_L005_R1_001_contig_2.g7). If that is a correct assumption, this should work for you:
$ FastaToTbl file.fa | sed 's/\.[^.]*\t/\t/' |
    awk -F"\t" '{
                    if(length($2)>length(a[$1])){
                        a[$1]=$0
                    }
                }
                END{
                    for(i in a){
                        print a[i]
                    }
                }' | TblToFasta
>Doug_NoIndex_L005_R1_001_contig_2.g7
atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgat
tcgtatgatcaactcatcattaatataaccaaaaacctagaaattctagccttcgatgat
gttgcagctgcggttcttgaagaagaaagtcggcgcaagaacaaagaagatagaccg

One possible issue with that approach is that it needs to keep one sequence per ID in memory. If you have a huge fasta file and not much memory, that might be a problem. If that's the case, you can first sort the file to ensure that similar IDs always appear together and then print each line as soon as you find the next ID:
FastaToTbl file.fa | sed 's/\.[^.]*\t/\t/' | LC_ALL=C sort -t$'\t' -k1,1 |
    awk -F"\t" '{
                    if(NR>1 && prev!=$1){
                        print a[$1];
                        prev=$1
                    }
                    if(length($2)>length(a[$1])){
                        a[$1]=$0
                    }
                }
                END{
                    print a[$1]
                }' | TblToFasta

FastaToTbl
#!/usr/bin/awk -f
{
        if (substr($1,1,1)==">")
        if (NR>1)
                    printf "\n%s\t", substr($0,2,length($0)-1)
        else
            printf "%s\t", substr($0,2,length($0)-1)
        else
                printf "%s", $0
}END{printf "\n"}

TblToFasta
#! /usr/bin/awk -f
{
  sequence=$NF

  ls = length(sequence)
  is = 1
  fld  = 1
  while (fld < NF)
  {
     if (fld == 1){printf ">"}
     printf "%s " , $fld

     if (fld == NF-1)
      {
        printf "\n"
      }
      fld = fld+1
  }
  while (is <= ls)
  {
    printf "%s\n", substr(sequence,is,60)
    is=is+60
  }
}


text:
Getting structure and sequence related to PDB IDs
I have two kinds of interactions: transient and stable. We are supposed to work on stable interaction, like interactions between two monomers in a heterodimer.
In a heterodimer there are two chains in which there are some residues between monomers which take part in the interaction of monomers and building the structure of heterodimer, and there are other residues interacting inside each monomer for building the structure of that monomer. There are some PDB IDs related to heterodimers in PDB.
Now imagine that we are trying to show which residues are physically interacting between two monomers in a heterodimer, and then we are trying two show the site of these physically interacting residues in the sequence related to that special structure. For that, we need to find the sequence related to that structure and we can download it from PDB.
But if we want to make a multiple sequence alignment of that sequence, first we need to find the refseq sequence of that and then we must blast the sequence and find homologues of that sequence.
Here there are some challenges:

we must find residues interacting in the structure of heterodimer between two monomers.
we must find the sequences of monomer chains of the heterodimer and map the interacting residues in the structure to the sequence.
we know that when a structure is solved and its monomers are sequenced, may be they are not able to completely sequence the heterodimer and residue numbers are maybe different from the related refseq sequence.

I would like to know how can we find the related refseq sequence (or sequences, in the case of heterodimer) and then how can we map the physically interacting residues in the structure of PDB ID to refseq sequence (finding the sites of these interacting residues in that refseq sequence) regarding three challenges mentioned above.
Imagine that I have 100 PDB IDs of 100 nonredundant heterodimers and I would like to find and download the structures, sequence and finally refseq sequence. What should I do?

So you want to map IDs to RefSeq IDs? Do you want protein or gene IDs?   Have you tried anything? If so, please [edit] and tell us what you tried and how it failed so we don't just give you the same solutions. What species are the PDBs from? What to PPIs have to do with it? PDB isn't a database of interactions. Could you [edit] your question and add a bit more detail so we can understand what you need, please?

You can find mappings PDB <-> Uniprot and Uniprot <-> RefSeq, but if you just want to use Blast, can't you use the sequence from PDB?

That works but I do not just want to do blast. I want to do multiple sequence alignment after that (after finding homologous sequences by running blast). The main challenge is mapping residue-residue interaction in PDB structure to refseq sequences and finding residues sites which participate in interaction. what can I do for that?

I think you could use the sequence from PDB also in MSA, and the alignment will show what the residues of interest correspond to in other sequences (but take what I write with a grain of salt, I don't have first-hand experience here)

You can download seqs and structures based on a list of PDB ids using http://www.rcsb.org/pdb/download/download.do#FASTA

The post might seem old, but I wish to share my thaughts based on how I solved similar problem. Here it goes:
I will use an example to illustrate. Other posts have already explained how you could get the protein ID from PDB and so on. I used the database PDBSUM (www.ebi.ac.uk/pdbsum/) and with protein PDB ID: 1KJF.
By following these steps I got to the results below:
Go to --> www.ebi.ac.uk/pdbsum/ -->(put in PDB ID) or (others e.g UniProt ID) --> Enter --> Click on "Ligands" -->(a LIGPLOT appears at the botoom of the page).
You can as well search just by pasting in the sequence.
Results:
You can see that 1KJF is made up of 3 chains A, B & P and both chains have clearly anntated and interacting residues. You can now see clearly read inidividual types of interactions.
Another good option is http://protein.bio.unipd.it/ring. It shows you a graph plot of all inteacting residues and how they are linked to each other. You can adjust the default result settings to meet your needs.
Hope it helps. :)

Yes.  'pdbsum' and 'ring' are definitely two useful tools for fast analytical work.