text:
Fast way to count number of reads and number of bases in a fastq file?
I am looking for a tool, preferably written in C or C++, that can quickly and efficiently count the number of reads and the number of bases in a compressed fastq file. I am currently doing this using zgrep and awk:
zgrep . foo.fastq.gz |
     awk 'NR%4==2{c++; l+=length($0)}
          END{
                print "Number of reads: "c;
                print "Number of bases in reads: "l
              }'

The files I need to process are regular ASCII text files (fastq) compressed using gzip (usually GNU gzip, possibly BSD gzip sometimes if clients are using macs).
zgrep . will print non-blank lines from the input file and the awk 'NR%4==2 will process every 4th line starting with the second (the sequence).
This works fine, but can take a very long time when dealing with large files such as WGS data. Is there a tool I can use (on Linux) that will give me these values? Or, if not, I'm also open to suggestions for speeding up the above command.
I know that the FASTQ specification (such as it is) allows for line breaks in the sequence and qual strings, so simply taking the second of every group of 4 lines is not guaranteed to work (see here). That's another reason why I'd rather use a more sophisticated, dedicated tool. However, I have yet to actually encounter a file with >4 lines per record in the wild, so I am not too worried about that.

It's difficult to get this to go massively quicker I think - as with this question working with large gzipped FASTQ files is mostly IO-bound. We could instead focus on making sure we are getting the right answer.
People deride them too often, but this is where a well-written parser is worth it's weight in gold. Heng Li gives us this FASTQ Parser in C.
I downloaded the example tarball and modified the example code (excuse my C...):
#include <zlib.h>
#include <stdio.h>
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

int main(int argc, char *argv[])
{
    gzFile fp;
    kseq_t *seq;
    int l;
    if (argc == 1) {
        fprintf(stderr, "Usage: %s <in.seq>\n", argv[0]);
        return 1;
    }
    fp = gzopen(argv[1], "r");
    seq = kseq_init(fp);
    int seqcount = 0;
    long seqlen = 0;
    while ((l = kseq_read(seq)) >= 0) {
        seqcount = seqcount + 1;
        seqlen = seqlen + (long)strlen(seq->seq.s);
    }
    kseq_destroy(seq);
    gzclose(fp);
    printf("Number of sequences: %d\n", seqcount);
    printf("Number of bases in sequences: %ld\n", seqlen);
    return 0;
}

Then make and kseq_test foo.fastq.gz.
For my example file (~35m reads of ~75bp) this took:
real    0m49.670s
user    0m49.364s
sys     0m0.304s

Compared with your example:
real    0m43.616s
user    1m35.060s
sys     0m5.240s

Konrad's solution (in my hands):
real    0m39.682s
user    1m11.900s
sys     0m5.112s

(By the way, just zcat-ing the data file to /dev/null):
real    0m38.736s
user    0m38.356s
sys     0m0.308s

So, I get pretty close in speed, but am likely to be more standards compliant. Also this solution gives you more flexibility with what you can do with the data.
And my horrible C can almost certainly be optimised.

Same test, with kseq.h from Github, as suggested in the comments:
My machine is under different load this morning, so I've retested. Wall clock times:
OP: 0m44.813s
Konrad: 0m40.061s
zcat > /dev/null: 0m34.508s
kseq.h (Github): 0m32.909s
So most recent version of kseq.h is faster than simply zcat-ing the file (consistently in my tests...).

Wow. So the C approach is actually marginally slower than the zgrep/awk I used and significantly slower than Konrad's trick of switching to wc.  That's surprising. I guess there must be some fancy parsing going on inside kseq_read.

The kseq.h from the tarball is outdated. The one from github is faster. That said, decompression is the bottleneck as you have observed.

@user172818 will that always be the case, irrespective of the size of the file? I would have thought that the decompression would be less of a bottleneck than the processing of each line when we're dealing with millions of lines. Is that wrong?

I think with respect to speed of core utility base solutions, people often overlook just how optimised these programs are (which of course are themselves C binaries), they have been around for a long time and are not going away.

Decompression will always be the bottleneck, it's the same situation we saw with https://bioinformatics.stackexchange.com/questions/361/what-is-the-fastest-way-to-calculate-the-number-of-unknown-nucleotides-in-fasta

The following is more than twice as fast; however, wc counts newline characters as well. We thus need to subtract the line count from the base count (using Bash):
fix_base_count() {
    local counts=($(cat))
    echo "${counts[0]} $((counts[1] - counts[0]))"
}

gunzip -c "$file" \
| awk 'NR % 4 == 2' \
| wc -cl \
| fix_base_count

All the caveats from Simon’s comment apply: this assumes the “simple” FASTQ format, where each record consists of exactly four lines. I think this is true for all files produced by Illumina sequencers and downstream tools.

What did you test this on? I used a 2.6G exome fastq file and this approach (without the bash function, so even faster) took  1m6.754s while my version using awk to count everything took 1m14.228s. What did you measure to see a two-fold increase in speed?

@terdon I’ve used the time command for benchmarking, and working on a 2.2GiB FASTQ file. The times on my machine were 9m5.628s and 3m5.825s real time, respectively. The user time in both cases was slightly elevated (11m vs 5m) to make the  difference a factor 2.

@terdon But given Simon’s answer and common sense, this is a fishy result indeed. I don’t know what my awk version is doing here.

pigz | awk | wc is the fastest method
First off for benchmarks with FASTQ it's best to use a specific real-world example with a known answer.  I've chosen this file:
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01815/sequence_read/ERR047740_1.filt.fastq.gz
as my test file, the correct answers being:
Number of reads: 67051220
Number of bases in reads: 6034609800

Next we want to find the fastest way possible to count these, all timings are the average wall-clock time (real) of 10 runs collected with the bash time on an otherwise unloaded system:
zgrep
zgrep . ERR047740_1.filt.fastq.gz |
     awk 'NR%4==2{c++; l+=length($0)}
          END{
                print "Number of reads: "c;
                print "Number of bases in reads: "l
              }'

This is the slowest method with an average run-time of 125.35 seconds
gzip awk
Using gzip we gain about another 10 seconds:
gzip -dc ERR047740_1.filt.fastq.gz |
     awk 'NR%4==2{c++; l+=length($0)}
          END{
                print "Number of reads: "c;
                print "Number of bases in reads: "l
              }'

Average run-time is 116.69 seconds
Konrad's gzip awk wc variant
fix_base_count() {
    local counts=($(cat))
    echo "${counts[0]} $((${counts[1]} - ${counts[0]}))"
}

gzip -dc ERR047740_1.filt.fastq.gz \
    | awk 'NR % 4 == 2' \
    | wc -cl \
    | fix_base_count

This runs slower on this test file than the gzip awk variant of the solution, average run-time is 122.28 seconds.
kseq_test using latest kseq.h from klib
Code compiled with: gcc -O2 -o kseq_test kseq_test.c -lz where kseq_test.c is Simon's adaptation of Heng Li's FASTQ parser.
kseq_test ERR047740_1.filt.fastq.gz
Average run-time is 99.14 seconds, which is better than the gzip core utilities based solution so far, but we can do better!
piz awk
Using Mark Adler's pigz as a drop-in replacement for gzip, note that pigz gives us a speed gain as on top of gzip as in addition to the main deflate thread it uses another 3 threads for reading, writing and checksum calculations, see the man page for details.
pigz -dc ERR047740_1.filt.fastq.gz |
     awk 'NR%4==2{c++; l+=length($0)}
          END{
                print "Number of reads: "c;
                print "Number of bases in reads: "l
              }'

Average run-time is now 93.86 seconds, this is ~5 seconds faster than the kseq based C code but we can further improve the benchmark.
pigz awk wc
Next we use pigz as a drop in replacment for Konrad's wc variant of the awk based solution.
fix_base_count() {
    local counts=($(cat))
    echo "${counts[0]} $((${counts[1]} - ${counts[0]}))"
}

gzip -dc ERR047740_1.filt.fastq.gz \
    | awk 'NR % 4 == 2' \
    | wc -cl \
    | fix_base_count

Average run-time is now down to 83.03 seconds, this is ~16 seconds faster than the kseq based solution and ~42 seconds faster than the OPs zgrep based solution.
Next as a baseline lets see just how much of this run-time is due to decompression of the input fastq.gz file.
gzip alone
gzip -dc ERR047740_1.filt.fastq.gz > /dev/null
Average run-time: 105.95 seconds, so the gzip based solutions (which also includes zcat and zgrep as these are provided by gzip) are never going to be faster than kseq_test.
pigz alone
pigz -dc ERR047740_1.filt.fastq.gz > /dev/null
Average run-time: 77.66 seconds, so quite clearly the additional three threads for read, write and checksum calculation offer a useful advantage.  What's more this speed-up is greater when leveraging the awk | wc based solution, it's not clear why, but I expect this is due to the extra write thread.
Interestingly average CPU usage across all threads is quite revealing for the various answers, I've collated these stats using GNU time /usr/bin/time --verbose
zgrep based solution 133% - must be more than one thread somehow
gzip | awk based solution 99% - all gzip based solutions run single-threaded at 99% CPU usage
pigz | awk 147%
gzip | awk | wc 99% as with gzip
pgiz | awk | wc 155%
kseq_test 99%
gzip > dev/null 99%
pigz > dev/null 155%
Whilst the main deflate thread in pigz will run at 100% CPU load the extra 3 don't quite fully occupy additional cores to 100% (as is evidenced by average CPU usage of ~150%) they do however clearly result in reduced run-time.
I'm using Ubuntu 16.04.2 LTS**, my gzip, zcat, zgrep versions are all gzip 1.6 and pigz is version 2.3.1.  gcc is version 5.4.0
** I think my patch level is actually 16.04.4 but I've not rebooted for 170 days :p

Nice and thorough! Thank you. I don't know why your tests showed pigz to be faster when it was slower in mine. I suspect it might be because I was testing using files stored on an NFS volume so my times could well have depended more on the network latency than anything else. I'll repeat your tests using a local drive and report back.

OK, confirmed. On your test file and on a local drive, pigz was indeed faster than gzip or zcat | grep.

wow excellent. thanks

I get fairly quick results with my fastx-length.pl script, with the added bonus of being able to handle multi-line FASTQ files and displaying additional read-length QC statistics:
time zcat albacored_all.fastq.gz | /bioinf/scripts/fastx-length.pl > /dev/null
Total sequences: 301135
Total length: 283.902419 Mb
Longest sequence: 5.601 kb
Shortest sequence: 6 b
Mean Length: 942 b
Median Length: 999 b
N50: 111835 sequences; L50: 1.103 kb
N90: 245243 sequences; L90: 608 b

real    0m8,802s
user    0m16,584s
sys 0m0,260s

Versus the script you have provided:
zcat albacored_all.fastq.gz | awk 'NR%4==2{c++; l+=length($0)}
          END{
                print "Number of reads: "c;
                print "Number of bases in reads: "l
              }'
Number of reads: 301135
Number of bases in reads: 283902419

real    0m8,382s
user    0m10,216s
sys 0m0,332s

Cat to /dev/null for comparison:
time zcat albacored_all.fastq.gz > /dev/null

real    0m7,877s
user    0m7,856s
sys 0m0,020s

I suspect that something using bioawk might be a bit faster (and similarly FASTQ-compliant).

bioawk is not that efficient, at least in my tests: https://bioinformatics.stackexchange.com/a/961/292

I hava implemented seqtk_counts using kseq.h from  klib
Just a few line of Codes:
#include <stdio.h>
#include <zlib.h>
#include "kseq.h"

KSEQ_INIT(gzFile, gzread)

int main(int argc, char *argv[]){

gzFile fp;
kseq_t *seq;

int l = 0;

int64_t total = 0;
int64_t lines = 0;

if (argc == 1) {
    fprintf(stderr, "Usage: %s <fastq> <sample>\n", argv[0]);

    return 1;
}

fp = strcmp(argv[1], "-")? gzopen(argv[1], "r") : gzdopen(fileno(stdin), "r");
seq = kseq_init(fp);

while ((l = kseq_read(seq)) >= 0){
    total += seq->seq.l;
    lines += 1;
}

printf("%s\t%lld\t%lld\n", argv[2] ,(long long)lines, (long long)total);
kseq_destroy(seq);
gzclose(fp);
return 0;

}
Compile it:
gcc  -O2  seqtk_counts.c  -o  seqtk_counts  -Iklib  -lz

Usage:
seqtk_counts foo.fasq.gz foo
or
cat foo.fasq.gz | seqtk_counts  - foo

Thanks for the edit! :) But isn't this basically the same program as shown in  sjcockell's answer?

Using pyGATB
(I use the same file as in https://bioinformatics.stackexchange.com/a/400/292, same workstation as in https://bioinformatics.stackexchange.com/a/380/292)
$ time python3 -c "from gatb import Bank; seq_lens = [len(seq) for seq in Bank('SRR077487_2.filt.fastq.gz')]; print('Number of reads: %d' % len(seq_lens), 'Number of bases in reads: %d' % sum(seq_lens), sep='\n')"
Number of reads: 23861612
Number of bases in reads: 2386161200

real    0m41.122s
user    0m40.788s
sys     0m0.312s

It is quite faster than bioawk:
$ time bioawk -c fastx '{nb_seq+=1; nb_char+=length($seq)} END {print "Number of reads: "nb_seq"\nNumber of bases in reads: "nb_char}' SRR077487_2.filt.fastq.gz
Number of reads: 23861612
Number of bases in reads: 2386161200

real    1m3.182s
user    1m2.916s
sys     0m0.268s

But not so much than the OP example:
$ time zgrep . SRR077487_2.filt.fastq.gz | awk 'NR%4==2{c++; l+=length($0)} END{print "Number of reads: "c; print "Number of bases in reads: "l}'
Number of reads: 23861612
Number of bases in reads: 2386161200

real    0m47.127s
user    1m36.292s
sys     0m6.796s

Or than the wc based solution:
$ fix_base_count() {
>     local counts=($(cat))
>     echo "${counts[0]} $((${counts[1]} - ${counts[0]}))"
> }
$ time gunzip -c SRR077487_2.filt.fastq.gz | awk 'NR % 4 == 2' | wc -cl | fix_base_count
23861612 2386161200

real    0m44.915s
user    1m12.000s
sys     0m6.972s

I didn't compare with C-based solutions.
The zcat to /dev/null reference is the following:
$ time zcat SRR077487_2.filt.fastq.gz > /dev/null

real    0m39.745s
user    0m39.464s
sys     0m0.252s

I'm still impressed by pyGATB speed

if the data is in SRA, there is sra-stat utility that returns reads,bases and quality distribution. these are stored in the SRA file. use --quick to get the stored stats or --statistics to calculate additional values broken down per readgroup/barcode.
sra-stat --quick SRR077487

Use assembly-stats!
$ assembly-stats barcode13_filtered.fastq
stats for barcode13_filtered.fastq
sum = 2080834976, n = 656192, ave = 3171.08, largest = 15321
N50 = 3598, n = 225863
N60 = 3263, n = 286569
N70 = 2920, n = 353904
N80 = 2548, n = 429982
N90 = 2081, n = 519671
N100 = 500, n = 656192
N_count = 0
Gaps = 0

real    0m7.418s
user    0m5.976s
sys     0m1.275s

Works on any fasta/fastq files, even those that aren't assemblies. And it's pretty fast!
Download from here: https://anaconda.org/bioconda/assembly-stats

Thanks, that looks interesting. How does the time compare to the existing answers here? Is this faster? And what are the numbers reported? Is that 656192 reads? And 2080834976 bases?

Yes, yes, and yes! 656192 reads, 2080834976 bases, with the respective N50 etc values. And compared to the ones above, it is lightning fast.
real 0m7.418s
user 0m5.976s

Could you show an example of that? You are running this on a single, decompressed file here, while the question is asking about compressed files and the other answers have shown that the rate limiting step is usually the decompression. So how fast is this approach when run on a compressed file and taking the decompression time into account? Also, in order for the reported times to be meaningful, please also include the time taken by one of the other methods (you can't compare times across different machines and files).

This should too good to be true, is it because the fastq is uncompressed? @sjcockell benchmark is on a comparable ~2.6e9 bases file, but it's about 5x slower... But it's a compressed file that was used.


text:
How can I count the number of reads that support a variant in a bam file?
I am calling variants from a human sample using bwa mem to align the reads and gatk to call the variants. I'm trying to understand why a specific variant was not called in my sample. I have checked the bam alignments in a GUI viewer and I can see that there are reads supporting the missing variant. It looks like the issue is a low allelic balance, with far more reads supporting the reference than the alternate allele but I want to get the actual numbers.
So, given a specific variant like this:
chr22   425236  C T

How can I count the number of reads in my sample.bam file that support that variant and the number that don't on Linux?

If you don't mind a bit of manual counting, then samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam will produce output where you can count the bases for that position. You could, of course, use the command line to do most of that automatically:
samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam |
    cut -f 5 | tr '[a-z]' '[A-Z]' | fold -w 1 | sort | uniq -c

That'll give you a count of how many of each base were seen.

ASCIIGenome (I'm the author) has a command, filterVariantReads, designed to inspect reads having a variant at a position or range. It would go along these lines:
ASCIIGenome -fa genome.fa aln.bam

Then go to the region of interest and use:
goto chr9:4917981-4918161
filterVariantReads -r 4918011

From this:

You get:

You can also print to screen or save to file the underlying sam records with:
print > reads.sam

For several regions, the whole process can be scripted for automation.
Hope this helps!

Hello, I am having a trouble to download this tool neither with conda nor from the source code. The .jar file is missing...

@user3224522 What have you tried? The zip file ASCIIGenome-1.16.0.zip should contain what you need. conda install asciigenome should also work. Please submit an issue to GitHub or post a question with more detail.

The variant record in the VCF produced by GATK should include information about depth -- both total read depth (DP) and allelic depth (AD), ie depth per allele.
Keep in mind however that the caller makes its decision based on more than just allelic depth. There are some additional statistics involved, including the quality of the individual bases. Also, if any filtering was applied, your variant may have been excluded due to other statistics, like for example if there is evidence of strand bias -- eg if all the reads supporting your variant are in the same orientation, that is a classic red flag suggesting it's an artifact.

source: Stack Exchange

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

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

sort: votes


text:
Is there a database of disease prevalence?
I want to collect data on how common different diseases are in different populations (or, at least, globally). Ideally, it would give me a way of querying the database with a disease name, and would return the number of cases per N population.
My main objective is to have a general idea of how prevalent a disease is, I am not trying to track current outbreaks, so I am not too bothered if the information is a little out of date, as long as the database is still maintained and not abandoned.
I have found FIDD, but it seems to be no longer maintained. Are there any actively maintained, similar tools that have the information I am after?

Comment as an answer, as requested.
The Incidence and Prevalance Database (IPD - subscription required).

For rare diseases, you can use Orphanet. You can download an xml file with prevalence from orphadata.
There is an OWL version of Orphanet called ORDO, you can browse it on the EBI Ontology Lookup Service. See for example Arrhythmogenic right ventricular cardiomyopathy and point prevalence relationships in the lower right box. There are a variety of programmatic ways to access this, including the OLS web API, SPARQL queries, or downloading the OWL and using an ontology API.

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

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


text:
Why does GATK produce both 0/1 and 1/0 genotypes in the same file? Are the two not equivalent?
I have always thought that 1/0 and 0/1 in VCF genotype fields are equivalent. And yet, GATK uses both. For example, these are two variants called in the same sample and the same run of GATK 4.1.4.0:
chr7    117120317   .   ATTCATTGTTTTGAAAGAAAGATGGAAGAATGAACTGAAG    A   748.97  .   AC=1;AF=0.5;AN=2;DP=64;ExcessHet=3.0103;FS=0;MLEAC=1;MLEAF=0.5;MQ=60;QD=11.89;SOR=7.223 GT:AD:DP:GQ:PL:SB   1/0:0,36:63:99:2294,1042,933:0,0,0,36
chr7    117120306   .   GACGCGCT    G   516.6   .   AC=1;AF=0.5;AN=2;BaseQRankSum=0.826;DP=70;ExcessHet=3.0103;FS=0;MLEAC=1;MLEAF=0.5;MQ=60;MQRankSum=0;QD=8.2;ReadPosRankSum=-5.786;SOR=0.976  GT:AD:DP:GQ:PL:SB   0/1:37,26:63:99:524,0,1203:0,20,0,26

And here are the same lines redacted for clarity:
chr7    117120317   .   AT..AG  A   [...]   GT:AD:DP:GQ:PL:SB   1/0:0,36:63:99:2294,1042,933:0,0,0,36
chr7    117120306   .   GA..CT  G   [...]   GT:AD:DP:GQ:PL:SB   0/1:37,26:63:99:524,0,1203:0,20,0,26

As you can see, one is genotyped as 1/0 and the other as 0/1. The VCF specs say:

GT (String): Genotype, encoded as allele values separated by either of / or |. The allele values are 0 for the
reference allele (what is in the REF field), 1 for the first allele listed in ALT, 2 for the second allele list in
ALT and so on. For diploid calls examples could be 0/1, 1 | 0, or 1/2, etc. Haploid calls, e.g. on Y, male
non-pseudoautosomal X, or mitochondrion, are indicated by having only one allele value. A triploid call might
look like 0/0/1. If a call cannot be made for a sample at a given locus, ‘.’ must be specified for each missing
allele in the GT field (for example ‘./.’ for a diploid genotype and ‘.’ for haploid genotype). The meanings of
the separators are as follows (see the PS field below for more details on incorporating phasing information into
the genotypes):

/ : genotype unphased
| : genotype phased

This suggests that although 1|0 is not the same as 0|1, the unphased 0/1 and 1/0 genotypes are completely equivalent. Which makes perfect sense given that the N/N genotypes are unphased, so n/N must be equivalent to N/n.
So, my questions are: why does gatk produce some lines with 0/1 and others with 1/0? Is it a bug? Am I missing something? I can't imagine a developer writing code to pick either 0/1 or 1/0 randomly, so presumably they indicate some difference in what gatk sees about a variant. So are 0/1 and 1/0 not 100% equivalent somehow?

 I have seen How is the GT field in a VCF file defined? which suggests the two are equivalent, but this seems to be contradicted by GATK's own output which is why I am asking this again to get a clear and preferably referenced answer.

@M__ this was from a human sample. But you can't know dominant or recessive unless you also know what gene, if any, is affected by the variant and what that gene's mode of inheritance is with respect to a specific disease, so I don't think it can be about AD vs AR.

Okay finally got it. The convention is it should be in ascending ordering, i.e. 0/1. The question is, why the convention broken?

I just realized that's almost certainly on me. I just took a couple of variants from a file randomly, and I hadn't really noticed the positions. I'm pretty sure they were in the right (ascending) order in the original file. Which lends more credence to your suggestion, @M__.

Convention
The convention is it should be in ascending ordering, i.e. 0/1. The question is, why the convention broken?
Phasing I think, at a guess, the answer is it depends what reference genome you are using. I would suggest that this result would not occur if the reference genome is mum (or dad) and the alleles are occurring in the child. The output supports this idea because it is marked with '/', i.e. they are not phased, thus the reference genome is neither mum nor dad.
Example provided The thing I would like to point out is in the example given 117120306 and 117120317 are not randomly scattered across the genome but appear to be 11 bp apart. For a human genome those are very strongly linked.
Linkage hypothesis Thus I would suggest the answer is linkage in this example, given the reference genome isn't a parent. The mutations are occurring on different strands therefore better represented as 0/1 and then 1/0, i.e. one strand comprises 'A' allele and the opposite strand is the 'G' allele.
Why is unphased linkage important In this scenario the interpretation could be hugely important because it would be stating that each allele has come from different parents. This would be more important than maintaining convention. In an inherited disorder there is a risk the child has a dominant phenotype, when each parent could have been heterozygous. Clearly it depends whether the mutation is forcing an amino acid mutation, at 11bp it is certainly possible, i.e. mutation in the first codon upstream and mutation in the second codon downstream.
More generally in any population analysis, linkage is an enormously important issue and central to any calculation. The point is simply that most analysis infers independent segregation of each allele and in this example that remains a good assumption because the allele (I am suggesting) are on opposite strands. If they were on the same strand it would create a bias away from panmixia.
The hypothesis is formed on this example alone and if the alleles were 10 megabp apart it would be impossible to conclude they were linked, therefore the convention of 0/1 should hold. That's my guess, its not a bad guess.

This is an interesting suggestion. I just picked two random variants from the file, so I will have to investigate a bit to see if what you suggest is actually the case more generally. I'd love to see this backed up by a gatk source (post, or dev, or source file or whatever). Certainly plausible, anyway, thanks!

Hope it helps @terdon. Its a hypothesis to work on. Geraldine@Broad will know of course.

I'm honored by your confidence in me, @M__ — But in this case I don't know for sure  I don't recall seeing something like this before. And if it was a phasing thing I would expect there to be phase related annotations. Are you seeing this in more than one record? And is this single or multi sample, have you done any subsetting, etc?
FWIW you might get a definitive answer faster on the GATK forum.

Alas, but regardless my confidence remains @Geraldine_VdAuwera  Thanks.

Thanks, @Geraldine_VdAuwera, I'll try there.


text:
How can I locate duplicated regions in a sequence?
I am facing an issue when trying to align short reads against a region in human chr5. The two Sensory Motor Neuron genes, (SMN1 and SMN2) are almost 100% identical and this causes the aligner to fail to align reads correctly since each read matches perfectly to two separate regions of the reference, leading to very low alignment scores.
If I BLAT the SMN1 gene's full genomic sequence against the human genome (hg19) on UCSC, I get (only showing the 1st two lines):
QUERY       SCORE START   END QSIZE IDENTITY  CHROM  STRAND  START       END   SPAN
--------------------------------------------------------------------------------------------------------------
NG_008691.1 35072     1 35072 35072   100.0%  chr5   +    70215768  70250839  35072
NG_008691.1 34279     1 35072 35072    99.8%  chr5   +    69340350  69375383  35034

And on the browser:

To get around this, I want to mask the SMN2 gene. However, the region of similarity extends beyond the annotated boundaries of the gene, so I would like to identify the full extent of this duplication.
So, given a large nucleotide sequence in fasta format, how can I compare it to itself and identify identical or near identical regions?

If you know approximately where are the duplicated regions, you can cut them out from the genome and align the two block using tools like MUMmer or Mauve. Although Mauve is not scalable solution (You can not scan genomes for duplicates with it) it's very handy for aligning small(ish) blocks. It can also align multiple sequences if you would like to check some tri/multiplications.
I downloaded the sequence of ch5 you linked and extracted big regions around the genes (300bp blocks). The first block is from 70100000 - 70400000 (with SMN1 sequence at 115768 - 150839) and the second from 69200000 - 69500000 (with SMN2 sequence at 140350 - 175383).
I used progressive mauve to align the two (click on File -> align with ProgressiveMauve -> add Sequence ... -> Align).

I found that the duplicated region is almost 200k long and maps between 44690 - 225520 of the first block and 19760 - 200940 of the second block, which we still need translate to the original genome coordinates, but I believe I can leave this to you :-).
P.S. you can also somehow upload annotations to Mauve, but I was lazy, so I "hand drawn" the arrow where I knew the genes are.

@KamilSJaron I can't get this to work. It keeps running out of memory even on my laptop which has 32G, of which 20 were available. I even tried aligning 2 files of only 100bp each and it crashed. I managed to get it to run on a server with 41G of RAM available, but I can't get the graphical output that way. I see you ran it on macOS, so it seems like a bug of the Linux version :(

GeneBlocks is a Python library for comparing DNA sequences. It can be used to find common blocks (identical regions) in a group of DNA sequences. Using a modified version of your example (I duplicated a small region):
import geneblocks

sequence = geneblocks.load_record('sequence.gb', name='seq')
sequence2 = geneblocks.load_record('sequence.gb', name='seq2')
blocks = geneblocks.CommonBlocks.from_sequences([sequence, sequence2])

blocks.plot_common_blocks()

Note that it wasn't designed to plot only one sequence, so as a workaround the sequence is given twice, with different names. As an alternative, you can use one sequence and plot with DNA Features Viewer:
import dna_features_viewer

blocks = geneblocks.CommonBlocks.from_sequences([sequence])
seq_to_plot = blocks.sequences_with_annotated_blocks()["seq"]
graphic_record = dna_features_viewer.BiopythonTranslator().translate_record(seq_to_plot)
ax, _ = graphic_record.plot(figure_width=10, strand_in_label_threshold=7)

See the documentation for more details. A web interface is also provided at EGF CUBA: find common blocks.
Disclaimer: I'm the current maintainer of GeneBlocks


text:
COG Annotation - Dealing with genes assigned to two or more COG categories
I am dealing with a list of genes which have been selected from a gene enrichment analysis. In order to see what kind of genes are overrepresented, I ran eggnog-mapper to do an orthology assignment of each gene against a bacterial database (provided in eggnog).
After running eggnog you get a .csv file containing for each gene which is the closest ortholog, and a list of GO, KEGG and COG terms associated.
From there, I would like to do a Fisher-test (also considering Bonferroni correction) to highlight COG categories overrepresented in our dataset compared to our reference genome.
Finally my question: How do you deal with genes assigned to multiple COG categories (same would apply for GO terms). You count +1 to each of the COG categories assigned? Do you just choose one category per gene?

You add 1 to the count of each COG category. You are looking for over represented COG categories so you must count them all. Many genes will be assigned to multiple categories. In fact, a majority of human genes are given multiple GO cellular function annotations. I have worked considerably in this field and one of the surprises I found is how common protein multifunctionality is. But I digress.
The point here is that you want to see whether any COG categories are over represented in your test set. Whether a gene is annotated to one or to a hundred categories is irrelevant. You are counting genes per COG, so you need to count each of the COGs assigned to a gene separately.


text:
Trouble using biomaRt to retrieve hgnc symbols from Ensembl transcript ids
I have a matrix of gene counts which I'm going to use as input for DESeq.  Right now, each gene is labeled by its Ensemble transcript ID, but I'd like to convert these to their HGNC symbols before I input them into DESeq for analysis. I'm attempting to do this conversion using biomaRt, however, this is my first time using the program, and I'm running into issues.  Here is my current code:
library("biomaRt")
mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL", dataset = "mmusculus_gene_ensembl")
transcript.ids <- rownames(txi.kallisto$counts)
hgnc_symbols <- getBM(attributes = "hgnc_symbol", filters = "ensembl_transcript_id", values = transcript.ids, mart = mart)

Printing hgnc_symbols shows that it isn't what I want:
>ghnc_symbols
[1] hgnc_symbol
<0 rows> (or 0-length row.names)

Lastly, here is what the input transcript IDs look like:
> head(rownames(txi.kallisto$counts))
[1] "ENSMUST00000178862.1" "ENSMUST00000178537.1" "ENSMUST00000196221.1"
[4] "ENSMUST00000179664.1" "ENSMUST00000177564.1" "ENSMUST00000179520.1"

I may be missing something here, but why would you be looking for HGNC (Human Genome Organisation  Gene Nomenclature Committee) names for mouse genes?

You need to specify the number without the version. Instead of "ENSMUST00000178862.1" just "ENSMUST00000178862":
You can do this with one more line:
g <- gsub("\\..*", "", rownames(txi.kallisto$counts))
(hgnc_symbols <- getBM(attributes = c("hgnc_symbol", "chromosome_name", "ensembl_transcript_id"), filters = "ensembl_transcript_id", values = g, mart = mart))
##  hgnc_symbol chromosome_name ensembl_transcript_id
##1          NA              14    ENSMUST00000177564
##2          NA               6    ENSMUST00000178537
##3          NA               6    ENSMUST00000178862
##4          NA              12    ENSMUST00000179520
##5          NA              14    ENSMUST00000179664
##6          NA              14    ENSMUST00000196221

With these transcripts it doesn't seem that there is an equivalent hgnc_symbol. But is recognized, otherwise they are not recognized:
(hgnc_symbols <- getBM(attributes = c("hgnc_symbol", "chromosome_name", "ensembl_transcript_id"), filters = "ensembl_transcript_id", values = rownames(txi.kallisto$counts), mart = mart))
## [1] hgnc_symbol           chromosome_name       ensembl_transcript_id
##<0 rows> (or 0-length row.names)

However you might be interested in doing something like this:
human = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
getLDS(attributes = "ensembl_transcript_id",
       filters = "ensembl_transcript_id", values = g,mart = mouse,
      attributesL = "hgnc_symbol", martL = human)
## The query to the BioMart webservice returned an invalid result: the number of
## columns in the result table does not equal the number of
## attributes in the query. Please report this to the mailing list.


text:
What is a simple command line tool for doing Needleman-Wunsch pair-wise alignment on the command line
I have two DNA strings:

GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC

and

AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG

I want a tool that allows me to do something like this on the command line:
$ aligner  GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG

and receive an ASCII visualization of a pairwise alignment. Something like this would work:
GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC---------
---------AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG

I prefer a tool that is packaged on bioconda. A tool that also does protein and RNA sequences is even more preferable.

Is this a toy example? If so, maybe update it to another one. These sequences are low complexity and the actual correct alignment output is unclear to me

The sequences are not important. Edits are welcome.

To get the alignment in your example, you have to use zero (or very small) gap cost at the ends of sequences. That is not the standard NW.

I'm not sure that's true? To quote from the EMBOSS package manual:
"In a Needleman-Wunsch global alignment, the entire length of each sequence is aligned. The sequences might be partially overlapping or one sequence might be aligned entirely internally to the other. There is no penalty for the hanging ends of the overlap. In bioinformatics, it is usually reasonable to assume that the sequences are incomplete and there should be no penalty for failing to align the missing bases."

You are looking for the needle program from the EMBOSS suite. Available in bioconda.
http://emboss.sourceforge.net/apps/release/6.6/emboss/apps/needle.html
To read sequences from the commandline, you need so specify the format as asis. To get the output on the screen, you'll need -stdout and to use the default alignment parameters (gap penalty 10, extend_penalty 0.5) you'll need -auto. Thus your query above would be:
$ needle -asequence asis:GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC -bsequence asis:AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG -auto -stdout
########################################
# Program: needle
# Rundate: Tue 27 Nov 2018 10:48:50
# Commandline: needle
#    -asequence asis:GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC
#    -bsequence asis:AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG
#    -auto
#    -stdout
# Align_format: srspair
# Report_file: stdout
########################################

#=======================================
#
# Aligned_sequences: 2
# 1: asis
# 2: asis
# Matrix: EDNAFULL
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 56
# Identity:      37/56 (66.1%)
# Similarity:    37/56 (66.1%)
# Gaps:          18/56 (32.1%)
# Score: 181.0
#
#
#=======================================

asis               1 GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC---     47
                              |||||||||||||||||||||||||||||||||||||.
asis               1 ---------AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAG     41

asis              47 ------     47

asis              42 AGGAGG     47

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

Install biopython with conda and use the following script:
#!/usr/bin/env python
import sys
from Bio import Align
aligner = Align.PairwiseAligner()
aligner.mode = "local"
alignments = aligner.align(sys.argv[1], sys.argv[2])
print(alignments[0])

The output is then:
.GGA-GGAGGGAG--AAGGAGGGAGGGAAGA-GGAGGGAG--AAGGAGGGAGGC
.|-|-||||||||--|||-|||-|||||-||-||||||||--|||-|||-|||.
AG-AAGGAGGGAGGGAAG-AGG-AGGGA-GAAGGAGGGAGGGAAG-AGG-AGG.

There are a LOT of alignments with that same score and you can tweak the details of the mode by changing aligner.mode.

yes, this should be pretty fast too since the alignment work is done in C (rather than python)

I'd just use something like t_coffee and write a little bash function to run it on the strings you give:
aligner(){
  tmp=$(mktemp);
  args=("$@")
  for ((i=0; i<$#; i++)); do
    printf '>%s\n%s\n' $i ${args[i]} >> "$tmp"
  done
  t_coffee "$tmp"  2>/dev/null | /bin/grep -Pv '^CLUSTAL|^\s*$'
  rm "$tmp"

}

If you add that to your shell's initialization file (e.g. ~/.bashrc if using bash on Linux, ~/.profile if using bash on macOS), you can run align2seqs with as many sequences as you want to align:
$ aligner 'GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC' 'AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG'
0               GG--AGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC
1               nAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGG-GAAGAGGAGG
                    *********      *******  *********  * ** *  *
$ aligner 'GAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC' 'GGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG' GAGGGAGGGAAGAGGAGGGAGAAGGATCCGAGGAAGAGGAGG AGAAGGAGGGAGGGAAGAGGCCCGAGATATATAAGGATCCGAGGAA
0               GA-GG----GAGAAGGAGGGAGGGAAGAGGAGGGAGAAGG--AGGGAGGC
1               GGAGG----GAGGGAAGAGGAGGGAGAAGGAGGGAG-GGA--AGAGGAGG
2               GA-GG----GAGGGAAGAGGAGGGAGAAGGATCCGA-GGA--AGAGGAGG
3               AGAAGGAGGGAGGGAAGAGGCCCGAGATATAT--AA-GGATCCGAG-GAA
                    *    ***      **   **     *       *    * *

Here's another option: https://github.com/noporpoise/seq-align; I've used it in the past as it is small and has no dependencies (just make and you're set).

This tool actually prompted my question! Do you know if anyone has packaged it?

Ah, no... I didn't notice the last part of your question.

Could you edit and include an example of how the OP could use this tool to do what they ask for? That will help future readers of this question.


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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

Yes, there bwa-mem was published as a preprint

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

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

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

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

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

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

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

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

# Constants
MEM_MAPQ_COEFFICIENT = 30.0
MEM_MAPQ_MAXIMUM = 60

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

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

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

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

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

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

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

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

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

    # Return the final mapping quality value.
    return mapping_quality

```


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

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

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

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

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

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

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

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

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

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

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

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

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

I tried to clarify my explanations.


text:
How can I calculate loss of heterozygosity (LOH) in NGS sequencing data?
I'm analyzing a tumor sample and a healthy sample from the same patient. I want to use sciClone to look at tumor clonality. One input is the genomic regions that need to be excluded due to LOH. I have both WES and WGS.
Is there currently a "gold standard" for which tool will return the most robust regions?
The sciClone paper mentions VarScan2, APOLLOH, HMMScan? Are these still state of the art?
Is it feasible to calculate LOH programmatically without the need of external software?

I think you will need external software. What you have ended up with?

Answer from @fabio-marroni, converted from comment:
WARNING! Self-promotion! We published a freely-available tool for detection of LOH (and also redeuction fo homozygosity). Here you find the link to the paper and therein will find the link to the software. We compared with other tools that you might also try to see which one is the best performing. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5466433/


text:
How can I use annotations to remove variants not relevant to cancer risk?
I currently have ~180 whole germlines and around 10M SNPs/indels. I would like to build a predictive model using Machine Learning (ML) techniques to predict cancer risk according to these germline variants. The thing is, most of these 10M variants are not relevant to cancer, so first I have to annotate them and remove non-relevant SNPs/indels.
Which tools/measurements do you think are important to use in order to filter out irrelevant variants and keep only, let's say, around 10,000 SNPs/indels?
I guess the most obvious is CADD score, so I could keep those SNPs/indels with PHRED >= 10.
Which one would you suggest taking into account that downstream analysis would be ML algorithms that take as features variants that really contribute to cancer?

While your question is specific to cancerous germline mutations, I'd suggest you look at the COSMIC database of somatic mutations to include in your analysis.
There are other factors to include in this kind of analysis you're suggesting, such as predictive deleterious effects (PolyPhen for example can perform such predictions).
If you have 10M variants/indels, for the same form of cancer, then look for common variants, or maybe the frequency of variants identified across exons of a gene.
What you've asked is highly broad and we'd need more specific details about your cohort to suggest more specific approaches.

I don't think it will be possible to do what you ask, right now with current knowledge. Selecting variants relevant to cancer risk is still an open problem and usually requires quite a lot of human intervention.
You can use different measures of population frequency to filter common variants with the assumption that frequent variants won't be pathogenic.
You can use multiple in-silico predictors of pathogenicity to try to remove "benign" or "tolerated" variants, but predictors are far from perfect and you have to take into account that many of them are essentially measuring the same, so when combining their outputs you may be biasing your analysis towards nucleotide conservation, etc...
You can use in-silico predictors of splicing effects. These are even farther from perfect and return very noisy data.
You can take into account where the variants lie in the gene. In many genes, variants in the last few nucleotides won't be pathogenic, but may be important in others. In some genes only variants in some domains have an effect, while not in others.
You have to take into account if only one or both alleles of certain genes are affected. For that you may need phasing information. And obviously the list of genes.
Some genes have an impact on cancer risk when they are "deactivated" (or broken somehow) (tumor supressors) while others have an effect when their function is improved (oncogenes).
And for many many genes you'll simply not know if they have any effect in cancer risk at all.
You can get data from different databases and annotate your genomes to remove many variants but in the end you won't have very high quality data to enter to your ML algorithms. Maybe this will work for your project, but take into account that this "first step" of data cleaning is not easy or obvious with current biological knowledge.


text:
How can I create my own GO association file (gaf)?
This question is based on a question on BioStars  posted >2 years ago by user jack.
It describes a very frequent problem of generating GO annotations for non-model organisms. While it is based on some specific format and single application (Ontologizer), it would be useful to have a general description of the pathway to getting to a GAF file.
Note, that the input format is lacking a bit of essential information, like how it was obtained. Therefore, it is har to assign evidence code. Therefore, lets assume that the assignments of GO terms were done automagically.

I want to do the Gene enrichment using Ontologizer without a
  predefined association file(it's not model organism).
I have parsed a file with two columns for that organism like this :
geneA  GO:0006950,GO:0005737
geneB  GO:0016020,GO:0005524,GO:0006468,GO:0005737,GO:0004674,GO:0006914,GO:0016021,GO:0015031
geneC  GO:0003779,GO:0006941,GO:0005524,GO:0003774,GO:0005516,GO:0005737,GO:0005863
geneD  GO:0005634,GO:0003677,GO:0030154,GO:0006350,GO:0006355,GO:0007275,GO:0030528

I have downloaded the .ob file from Gene ontology file which contain
  this information (from here) :
!
! GO IDs (primary only) and name text strings
! GO:0000000 [tab] text string [tab] F|P|C
! where F = molecular function, P = biological process, C = cellular component
!
GO:0000001  mitochondrion inheritance   P
GO:0000002  mitochondrial genome maintenance    P
GO:0000003  reproduction    P
GO:0000005  ribosomal chaperone activity    F
GO:0000006  high affinity zinc uptake transmembrane transporter activity    F
GO:0000007  low-affinity zinc ion transmembrane transporter activity    F
GO:0000008  thioredoxin F
GO:0000009  alpha-1,6-mannosyltransferase activity  F
GO:0000010  trans-hexaprenyltranstransferase activity   F
GO:0000011  vacuole inheritance P

What I need as output is .gaf file in the following format (in the
  format of the files here):
!gaf-version: 2.0

!Project_name: Leishmania major GeneDB

!URL: http://www.genedb.org/leish

!Contact Email: mb4@sanger.ac.uk

 GeneDB_Lmajor    LmjF.36.4770    LmjF.36.4770        GO:0003723    PMID:22396527    ISO    GeneDB:Tb927.10.10130    F    mitochondrial RNA binding complex 1 subunit, putative    LmjF36.4770    gene    taxon:347515    20120910    GeneDB_Lmajor
 GeneDB_Lmajor    LmjF.36.4770    LmjF.36.4770        GO:0044429    PMID:20660476    ISS        C    mitochondrial RNA binding complex 1 subunit, putative    LmjF36.4770    gene    taxon:347515    20100803 GeneDB_Lmajor             GeneDB_Lmajor    LmjF.36.4770    LmjF.36.4770        GO:0016554    PMID:22396527    ISO    GeneDB:Tb927.10.10130    P    mitochondrial RNA binding complex 1 subunit, putative    LmjF36.4770    gene   taxon:347515    20120910    GeneDB_Lmajor
 GeneDB_Lmajor    LmjF.36.4770    LmjF.36.4770        GO:0048255    PMID:22396527    ISO    GeneDB:Tb927.10.10130    P    mitochondrial RNA binding complex 1 subunit, putative    LmjF36.4770    gene    taxon:347515    20120910    GeneDB_Lmajor

How to create your own GO association file (gaf)?

Please add some tags for GO, non-model organism, gaf, overrepresentation analysis etc. I just used a random tag because there are not many sensible tags in the list yet.

I modified a bit the question, added tags and I centered the question in creating the file in that format, if you disagree, just revert my changes

I added [tag:gene-ontology] and [tag:gaf]. I'm not sure whether a [tag:over-representation-analysis] or a [tag:enrichment-analysis] would be more appropriate.

@Llopis: Hopefully I didn't just accidentally overwrite your changes!

No luck, my edit was rejected (or at least doesn't appear) :\

If you want to bring questions from another site, please don't just copy/paste them verbatim. Instead, ask them again in your own words. And if you do copy, it would be better to put the quoted text in a quote block to clearly indicate that it is a quote.

@terdon, sorry need to adapt to the new rules ;) I accepted the edits, however I don't have time to rewrite the question in my own words, and I don't think it is strictly necessary. Even if that might seem a contradiction with my previous post, it is not imo, because the question is pretty aged.

@Michael no, I didn't expect you to rewrite this now! But please do so next time. Going around just pinching questions verbatim from other sites doesn't seem right. Not unless you're in a position to offer clarifications or this is an issue that affects you too.

Indeed, I can offer ample clarification for this question, we are doing exactly the same thing, except that we are using InterproScan to assign GO terms not Blast2GO as was likely used here. I have a similar script to convert IPS output to GAF.

Ha! I only now realized you were the person who first raised the whole copying from other sites issue as a problem :)

Here's a Perl script that can do this:
#!/usr/bin/env perl
use strict;
use warnings;

## Change this to whatever taxon you are working with
my $taxon = 'taxon:1000';
chomp(my $date = `date +%Y%M%d`);

my (%aspect, %gos);
## Read the GO.terms_and_ids file to get the aspect (sub ontology)
## of each GO term.
open(my $fh, $ARGV[0]) or die "Need a GO.terms_and_ids file as 1st arg: $!\n";
while (<$fh>) {
    next if /^!/;
    chomp;
    my @fields = split(/\t/);
    ## $aspect{GO:0000001} = 'P'
    $aspect{$fields[0]} = $fields[2];
}
close($fh);

## Read the list of gene annotations
open($fh, $ARGV[1]) or die "Need a list of gene annotattions as 2nd arg: $!\n";
while (<$fh>) {
    chomp;
    my ($gene, @terms) = split(/[\s,]+/);
    ## $gos{geneA} = (go1, go2 ... goN)
    $gos{$gene} = [ @terms ];
}
close($fh);

foreach my $gene (keys(%gos)) {
    foreach my $term (@{$gos{$gene}}) {
        ## Warn and skip if there is no aspect for this term
        if (!$aspect{$term}) {
            print STDERR "Unknown GO term ($term) for gene $gene\n";
            next;
        }
        ## Build a pseudo GAF line
        my @out = ('DB', $gene, $gene, ' ', $term, 'PMID:foo', 'TAS', ' ', $aspect{$term},
                             $gene, ' ', 'protein', $taxon, $date, 'DB', ' ', ' ');
        print join("\t", @out). "\n";
    }
}

Make it executable and run it with the GO.terms_and_ids file as the 1st argument and the list of gene annotations as the second. Using the current GO.terms_and_ids and the example annotations in the question, I get:
$ foo.pl GO.terms_and_ids file.gos
DB  geneD   geneD       GO:0005634  PMID:foo    TAS     C   geneD       protein taxon:1000  20170308    DB
DB  geneD   geneD       GO:0003677  PMID:foo    TAS     F   geneD       protein taxon:1000  20170308    DB
DB  geneD   geneD       GO:0030154  PMID:foo    TAS     P   geneD       protein taxon:1000  20170308    DB
Unknown GO term (GO:0006350) for gene geneD
DB  geneD   geneD       GO:0006355  PMID:foo    TAS     P   geneD       protein taxon:1000  20170308    DB
DB  geneD   geneD       GO:0007275  PMID:foo    TAS     P   geneD       protein taxon:1000  20170308    DB
DB  geneD   geneD       GO:0030528  PMID:foo    TAS     F   geneD       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0016020  PMID:foo    TAS     C   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0005524  PMID:foo    TAS     F   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0006468  PMID:foo    TAS     P   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0005737  PMID:foo    TAS     C   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0004674  PMID:foo    TAS     F   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0006914  PMID:foo    TAS     P   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0016021  PMID:foo    TAS     C   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0015031  PMID:foo    TAS     P   geneB       protein taxon:1000  20170308    DB
DB  geneA   geneA       GO:0006950  PMID:foo    TAS     P   geneA       protein taxon:1000  20170308    DB
DB  geneA   geneA       GO:0005737  PMID:foo    TAS     C   geneA       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0003779  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0006941  PMID:foo    TAS     P   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0005524  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0003774  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0005516  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0005737  PMID:foo    TAS     C   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0005863  PMID:foo    TAS     C   geneC       protein taxon:1000  20170308    DB

Note that this is very much a pseudo-GAF file since most of the fields apart from the gene name, GO term and sub-ontology are fake. It should still work for what you need, however.

I think the evidence code should be changed to IEA, because it is unlikely one would have manually annotated GO's for a non-model organism.

@Michael, it doesn't really make any difference so I chose TAS since these are annotations provided by the OP. But sure, IEA would also work. The field will be ignored in any case when doing an enrichment analysis so you can put whatever you prefer.

Sure, IEA might be even discarded by some tools.


text:
Verify a predicted protein in one genome in a different genome of the same species
I have two genome assemblies of the same non-model species, call them Assembly 1 (generated from Illumina data) and Assembly 2 (generated from PacBio data).
For Assembly 1, I also have predicted proteome data, generated with EVM. Say there is a protein, call it Protein X, for which I have the fasta sequence from predictions, but suspect an error in Assembly 1 (because the Protein X has a long segment of auto-identity). For this reason, I would like to verify the assembly of the locus coding for Protein X in Assembly 2 and/or verify the presence of Protein X in the Assembly 2.
Data in my possession:

both genome assemblies (fasta)
raw rna-seq data
reference-based transcriptome assembly for each genome assembly (fasta/gtf)
set of predicted protein sequence based on Assembly 1 (fasta)

Things I was thinking of:

build the EVM predictions on the Assembly 2 and verify the identity of Protein X by simple fasta identity check
look for the protein sequence in the Assembly 2. I could use something similar to Reverse Translate to get the most likely nucleotide sequence from the input amino-acid sequence of Protein X, then align it to the Assembly 2 to check the coordinates of the locus it originates from - but I am not very sure of the relevance and/or correctness of this procedure
something allowing 'blast' of the protein sequence in the Assembly 2 (but I am not aware of any software that could do that)

Any ideas on how I could do that?

This should be very easy to do. Here are some options:

Use a tool like Exonerate or GeneWise that can align protein sequences to genomic DNA while attempting to model splice sites etc.
As you said, blast. There's nothing special in what you describe, really, this is what tBLASTn is for. Just use your protein as the query sequence and tBLASTn as the blast flavor which will take a protein query and match it against all 6 possible reading frames of a DNA database. Set up your assembly as the blast DB, and then blast as you would normally.

Of the two, using the first approach is much better since it will try to build a valid gene model for you instead of simply finding the regions of high homology as blast will do.
Finally, just for completeness, you can go old school and run a full de-novo gene prediction tool on your assembly. Something like GeneID for instance. But that's certainly overkill if you're just looking for a single gene.


text:
What is this 5-column sequencing format from 2009?
There is high throughput sequencing data here, and I don't know what format it is in.
It was submitted in 2009, and the description says the following:

Library strategy:     ncRNA-Seq
Library source:   transcriptomic
Library selection:    size fractionation
Instrument model:     Illumina Genome Analyzer II
Description:  CIPPNK, tar file of Illumina *_seq.txt files provided as supplementary file

I got the archive here:
ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM455nnn/GSM455387/suppl/GSM455387%5FWT%5FCIPPNK%5Fseq%5Fs1%2Etar%2Egz
Inside, there are 330 files from s_1_0001_seq.txt to s_1_0330_seq.txt which are tab-delimited text files where the first column is always 1, the second has the number found in the file name, then 2 mysterious integers, and then what looks like a read of length 36, with sometimes a dot instead of a letter:
$ head s_1_0330_seq.txt
1   330 690 785 TTCCTACATTGTTCCCCCATGCTGTTGGCACCATCA
1   330 44  145 TTTTTATCACGAGTTTTAAATCTGTAGTCACCATCA
1   330 53  141 AATAATGCATAACAAAACGGAATCTGTAGAA.AAA.
1   330 784 461 TAATTGTAGTGATTGATCAATCTGTAGGCACCATCA
1   330 588 634 TATTATGCACATTTTCTAGTTCACTGTAGGCACCAT
1   330 718 678 TTACATGTTTCGGGTAGGAGCCTGTAGGCACCATCA
1   330 635 834 TGTGATCATTAGTTCAAAGCCCCCTGTCGGCACCCT
1   330 494 523 TGAAAATCAAAAATGCTGAACTGTAGGCACCATCAA
1   330 393 783 TTTTTTTTTAAATTTAAAAAAACTGTAGGCACCATC
1   330 48  148 GTTTAACCGTGTAGACGTTGGTTTCTGTAGGCACCA

This seems to be some "old" high-throughput sequencing format. I think someone in a 2008 message in sequanswers was dealing with this type of file:
http://seqanswers.com/forums/showpost.php?p=1841&postcount=8
What is this format that seemed so standard back then, that the authors did not give more info than describing the files as "Illumina *_seq.txt files"? I do not dare ask them for such a trivial question (the given contact is a Nobel prize laureate and is probably too busy to answer random bioinformatics questions).
In particular, what are the columns 3 and 4, and what do the dots mean?

This is an early Solexa/Illumina sequencer format. The columns are the identifier to location on the flowcell.  I believe the "." was the original placeholder for an unread base, which has been replaced by and "N" in current Illumina sequencing output.
From on http://www.crg.eu/en/content/processing-and-analysis-illumina-sequencing-data.

seq.txt (Gerald)
4 1 23 1566 ACCGCTCTCGTGCTCGTCGCTGCGTTGAGGCTTGCG
The first 4 columns are an ID that provides location details of the cluster in the flow cell, followed by the read sequence. Seq.txt file size typically ranges between 250 Megabytes and 1 Gigabyte. One file is generated per lane.

Ah, so if the authors did not provide other files, I guess the only way I can convert this to fastq is to chose arbitrary qualities.

Correct this format did not contain quality scores.  There was a qseq.txt file that did.  Although to be honest back then the quality scores were no that meaningful as the Illumina algorithm was not very accurate in assigning them.


text:
Parsing gtf file for transcript ID and transcript name
I used this previously to get ensembl ID and its gene name from the gencode gtf annotation file
This was the script I used
awk '{
    for (i = 1; i <= NF; i++) {
        if ($i ~ /gene_id|gene_name/) {
        printf "%s ", $(i+1)
        }
    }
    print ""
}' Homo_sapiens.GRCh37.70.gtf | sed -e 's/"//g' -e 's/;//g' -e 's/ /\t/' | sort -k1,1 | uniq > Homo_sapiens.GRCh37.70.txt

And the input file:
##description: evidence-based annotation of the human genome (GRCh38), version 21 (Ensembl 77)
##provider: GENCODE
##contact: gencode@sanger.ac.uk
##format: gtf
##date: 2014-09-29
chr1    HAVANA  gene    11869   14409   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENSG00000223972.5"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1"; level 2; havana_gene "OTTHUMG00000000961.2";
chr1    HAVANA  transcript  11869   14409   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_status "KNOWN"; transcript_name "DDX11L1-002"; level 2; tag "basic"; transcript_support_level "1"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    11869   12227   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_status "KNOWN"; transcript_name "DDX11L1-002"; exon_number 1; exon_id "ENSE00002234944.1"; level 2; tag "basic"; transcript_support_level "1"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    12613   12721   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_status "KNOWN"; transcript_name "DDX11L1-002"; exon_number 2; exon_id "ENSE00003582793.1"; level 2; tag "basic"; transcript_support_level "1"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    13221   14409   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_status "KNOWN"; transcript_name "DDX11L1-002"; exon_number 3; exon_id "ENSE00002312635.1"; level 2; tag "basic"; transcript_support_level "1"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  transcript  12010   13670   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000450305.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1-001"; level 2; ont "PGO:0000005"; ont "PGO:0000019"; tag "basic"; transcript_support_level "NA"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000002844.2";
chr1    HAVANA  exon    12010   12057   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000450305.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1-001"; exon_number 1; exon_id "ENSE00001948541.1"; level 2; ont "PGO:0000005"; ont "PGO:0000019"; tag "basic"; transcript_support_level "NA"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000002844.2";
chr1    HAVANA  exon    12179   12227   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000450305.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1-001"; exon_number 2; exon_id "ENSE00001671638.2"; level 2; ont "PGO:0000005"; ont "PGO:0000019"; tag "basic"; transcript_support_level "NA"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000002844.2";
chr1    HAVANA  exon    12613   12697   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000450305.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1-001"; exon_number 3; exon_id "ENSE00001758273.2"; level 2; ont "PGO:0000005"; ont "PGO:0000019"; tag "basic"; transcript_support_level "NA"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000002844.2";
chr1    HAVANA  exon    12975   13052   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000450305.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1-001"; exon_number 4; exon_id "ENSE00001799933.2"; level 2; ont "PGO:0000005"; ont "PGO:0000019"; tag "basic"; transcript_support_level "NA"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000002844.2";

What I'm trying to do now is replace the gene_ID with transcript_ID and the gene_name with transcript_name but I'm doing something wrong as I'm not getting the transcript_ID but the gene Ensembl ID instead. I'm not sure if I'm choosing the right column for parsing the transcript_ID and name.
So I'm expecting to get this
transcript_id "ENST00000450305.2";  gene_name "DDX11L1";transcript_name "DDX11L1-001"

Any help or suggestion would be highly appreciated

I get the Ensembl transcript ID and the transcript name (e.g., "ENST00000320153 DDI2-201"). What is the output you're getting?

im getting ENSG and ENST  ids both along with gene name, did you run the same thing?im doing something wrong perhaps

I'm running the exact same thing with /transcript_id|transcript_name/ rather than /gene_id|gene_name/.

okay ,yes im doing the same thing...but my output is everything combined gene as well as transcript ID with names if you can look into my file may be you can tell https://drive.google.com/open?id=1Frzy0R1j1Vncog6LP5E8N2TIJFBZsS3z what i need is can i get only transcript ID and name...

I don't see how your awk command would work since you're using whitespace as a field delimiter. In any case, you can use a much shorter perl command for this:
$ perl -lne 'print "@m" if @m=(/((?:transcript_id|gene_name|transcript_name)\s+\S+)/g);' Homo_sapiens.GRCh37.70.gtf
transcript_id "ENSG00000223972.5"; gene_name "DDX11L1"; transcript_name "DDX11L1";
transcript_id "ENST00000456328.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-002";
transcript_id "ENST00000456328.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-002";
transcript_id "ENST00000456328.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-002";
transcript_id "ENST00000456328.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-002";
transcript_id "ENST00000450305.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-001";
transcript_id "ENST00000450305.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-001";
transcript_id "ENST00000450305.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-001";
transcript_id "ENST00000450305.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-001";
transcript_id "ENST00000450305.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-001";

The scriptlet will print the contents of the @m array if this line matches what you need. The regex magic is:

(?: ): this is a non-capturing group. We need this so that @m only holds what is inside the capturing parentheses (see below).
(  (?:transcript_id|gene_name|transcript_name)\s+\S+) : the outer parentheses will capture everything matched: one of transcript_id|gene_name|transcript_name followed by 1 or more whitespace characters (\s+) and one or more non-whitespace characters.
print @m if @m=... : only print if the match was successful.

i was in hurry and dire need almost all the work i do it in R ,i tried it in R due to the size of the gtf it just froze ,  looks like I need to learn some basics of awk along with perl that will make my life simpler...thank you very much i will use your code...

that was super fast in takes time with awk but this is the time it took
real 0m16.112s
user 0m14.653s
sys 0m0.693s


text:
Extracting features or gene from PCA after calculating PCA for downstream analysis
The quote below is from this paper:

We performed principal component analysis (PCA) of low-coverage sequencing data to identify genes explaining variation across cells. PCA separated the cells into groups corresponding to the source populations (Fig. 2c and Supplementary Figs. 3–5), and the genes distinguishing each group reflected the biological properties of the cell types (Supplementary Fig. 5 and Supplementary Table 3). PCA of low- and high-coverage sequencing data revealed a remarkably similar graphical distribution of analyzed cells, and the majority (78%) of the top 500 genes determined by PCA were shared between PCA performed on low- and high-coverage data (Supplementary Figs. 4 and 6 and Supplementary Table 4).

How are they performing PCA and then finding out the components which contribute to the strongest PC and then taking out the genes for downstream analysis? Let's say there are genes involved in PC1 and PC2 which define a certain cell type and that distinguish it. I think for a biologist, like me, I just want to get the genes or gene list that are getting involved in determining lineage or cell type. How can I do this?
As a to test what they might be doing i did this
test<- read.csv('NON_CODING.csv',header = T,row.names = 1)

# preform PCA
pca = prcomp(t(test), center=TRUE, scale=TRUE)

then after calculating PCA i do see this in my console when I type pca
> pca
Standard deviations (1, .., p=16):
 [1] 3.826851e+01 2.080405e+01 1.568739e+01 1.349256e+01 1.119348e+01 9.980547e+00 8.365034e+00
 [8] 8.098841e+00 7.519507e+00 6.184505e+00 5.880260e+00 5.139609e+00 4.851091e+00 4.335606e+00
[15] 3.870918e+00 3.173450e-14

Rotation (n x k) = (2899 x 16):
                                       PC1           PC2           PC3           PC4
5S_rRNA                      -1.090574e-02 -2.412665e-03  1.637689e-02 -3.603865e-02
AB019441.29                  -1.928250e-02  1.821083e-03 -9.713724e-03 -1.978737e-03
ABBA01017803.1               -1.823266e-02 -3.727144e-03 -9.790131e-04 -3.937062e-02
ABC14-1080714F14.1            2.438024e-02  3.816019e-03  1.657784e-04 -7.141480e-03
ABC7-481722F1.1               2.467403e-02  2.873432e-03  3.781153e-03  3.790824e-03
AC000036.4                   -1.066777e-02  2.838305e-02  2.642673e-02 -7.998608e-03
AC000089.3                   -1.249026e-02 -7.870864e-04 -2.569602e-02 -3.165073e-03
AC000120.7                   -1.313696e-02  5.923245e-03  1.633255e-02  4.077602e-02
AC000123.4                    2.415996e-02  1.029732e-02  8.930792e-03 -1.118660e-02
AC000403.4                    8.692835e-03 -3.060170e-02  3.746839e-02 -6.968971e-03
AC002064.4                   -1.799590e-02  2.147160e-02 -2.303742e-02  7.647937e-03

Now how do i find which are the list or set of genes that gives me PC1 or PC2 so on
Any suggestion or help how to get the genes that make most difference between two major component that would be very helpful

The help page explains that there is a component which:

rotation
  the matrix of variable loadings (i.e., a matrix whose columns contain the eigenvectors). The function princomp returns this in the element loadings.

So using using the rotation one can use it to calculate the importance of each feature in the input:
C <- chol(S <- toeplitz(.9 ^ (0:31))) # Cov.matrix and its root
set.seed(17)
X <- matrix(rnorm(32000), 1000, 32)
Z <- X %*% C  ## ==>  cov(Z) ~=  C'C = S
pZ <- prcomp(Z)
genesImportance <- Z %*% pZ$rotation
genesImportance2 <- pZ$x
all.equal(genesImportance, genesImportance2)
## [1] "Mean relative difference: 0.03680457"

This new matrix has the loading for each gene for each component. So you can now select them according to how much do they participate in each component or whatever you want. Or you can use it to predict the position of samples you didn't set in the initial data.

Im using prcomp function would that work as you mentioned princomp ?

Z%*%pZ$rotation should be the same as pZ$x (that’s in fact how pZ$x is calculated by prcomp).

Totally true, I didn't went further in the help page, but with rotation one can predict the position of new samples not included in the original data. I'll edit it.

@KonradRudolph I'm still not getting where is this coming from  italic pZ$x   im little ignorant as I see examples and do it so i have to do a back referencing  to in R methods and tutorials , is it in my data frame as well ?with context to my data frame are you suggesting this  italic"pca$x"  and what is that  Z  actually

@krushnachChandra Z is the initial data for the PCA. pZ is the output of prcomp. See the edit, I might not explain it clearly enough

okay i will try to find it out

@Llopis i tried that i get error like this Error in Z %*% pZ$rotation :
requires numeric/complex matrix/vector arguments

@krushnachChandra See the new edit. Could you explain what code are you using and what is Z?