# Question
Title: 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.fasq.gz |
     awk 'NR%4==2{c++; l+=length($0)}
          END{
                print "Number of reads: "c;
                print "Number of bases in reads: "l
              }'

```

The `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.

# Answer
> 14 votes
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...).

# Answer
> 7 votes
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.

# Answer
> 7 votes
# 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

# Answer
> 5 votes
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).

# Answer
> 2 votes
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

```

# Answer
> 1 votes
### 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

# Answer
> 1 votes
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

# Answer
> 1 votes
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

---
Tags: ngs, fastq, software-recommendation, benchmarking
---


# Question
Title: Get frequency of index sequences of fastq file
I recently did some sequencing on a MiniSeq and unexpectedly the `Undetermined_S0_L001_R1_001.fastq.gz` file contained quite a lot of data (about 40% of total). I would like to know which indices are present in this file and ideally get a frequency of all index sequences in this file. How can I do this?

This is the output of `head`:

```
@MN00148:128:000H2YCW2:1:11102:12116:1065 1:N:0:0
GNAGAGGATTCAGGGAGCCAGCGCGCATATATCAG
+
F#FFFFFFFFFFFFFFFFFFA=/A///////////
@MN00148:128:000H2YCW2:1:11102:11354:1065 1:N:0:0
GNACCCGTGCGGCTGGAGAGACCACTAGAAAGAAT
+
A#AFFFFFFFFFFFFFFFFF///////FF=////6
@MN00148:128:000H2YCW2:1:11102:4103:1066 1:N:0:0
ANGTTGGTGGTACATCAGGGACGCGGAGATCTCAG

```

# Answer
> 1 votes
I actually went back to the output folder on the MiniSeq and found this file: `DemultiplexSummaryF1L1`, which contains demultiplexing information of the whole run. Half-way through this file I found what I was looking for:

```
### Most Popular Index Sequences

### Columns: Sequence ReverseComplement HitCount

Index
TGACCA  TGGTCA  2455964
CGATGT  ACATCG  2320579
GGCCCC  GGGGCC  2023652
CGCGGG  CCCGCG  1983605
CCGTCC  GGACGG  1641274
CAGATC  GATCTG  1514252
CCGGCC  GGCCGG  1319383
CCGCGC  GCGCGG  1109940
ACAGTG  CACTGT  891344
CCCGGG  CCCGGG  706083
TGCCCA  TGGGCA  168294
CAGCTC  GAGCTG  164091
.... [data continues]

```

# Answer
> 1 votes
I don't know that you can get them from that file. I've noticed that when the Miseq is left to run its version of fastq generation software, there are no indices in the undetermined file. I usually run bcl2fastq myself if I want that information.

They recently came out with a new version of software for the MiSeq, so things might be different now.

---
Tags: ngs, sequence-analysis, fastq, demultiplexing, indexing
---


# Question
Title: How to estimate the phenotypic variation explained by top SNPs from a GWAS study?
I have conducted a large-scale GWAS study and got a few significantly associated SNPs. I used `GEMMA` with `-lmm 1` options to run the GWAS and obtain the `beta` and `standard-error` estimates. I want to estimate the percent phenotypic variation explained by each of the significant SNPs. I used the following procedure for estimating the variance explained in `R`:

```
fit <- lm (Phenotypic_value ~ SNP_data, data = a)
summary(fit)$adj.r.squared

```

Here, the datafile `a` contains three columns namely, `sample_ID`, `Phenotypic_value` for each sample, and the biallelic `SNP_data`. I got a value which is 0.43 meaning 43% phenotypic variation explained by the SNP.

Again, I used another formula which is: `2*f*(1-f)*b.alt^2`. Here, `f` is the minor allele frequency and `b.alt` is the effect size i.e. `beta estimate` obtained from `GEMMA`. This gives me a value of 0.03 meaning 3% variation explained which seems reasonable to me.

My question is that which of the above method is correct? Why do I obtain such a huge difference between the two approaches? Is there any other way to estimate the percent variation explained?

Alternatively, from the GEMMA google group, I have got this formula `pve <- var(x) * (beta^2 + se^2)/var(y)`. But I do not understand how can I obtain the value of `var(x)` and `var(y)`.

I borrowed the first and the last formula from GEMMA google group discussion. The first one is obviously a normal linear model equation in statistics. The second formula is from link. Also, the data structure is three columns, 1 is with individual, 2 is the genetic information (individual carrying A or T) and 3 is the phenotype. In the last formula, `var(x)` is the variance of the genotype vector as they said in GEMMA group and the `beta` is the effect size estimate and `se` is the standard error from GEMMA output.

It will be great to receive some feedback on this. Thank you.

# Answer
*Answer from @maximilian-press, converted from comment:*

At a guess, the GEMMA linear mixed model is controlling for population structure and other SNPs, which your R example is not (it only fits 1 SNP?), but it's hard to be sure with this information. It would be helpful to see data samples. There is probably more context and usage information available in those resources beyond just the equations. It is especially difficult to evaluate with no terms defined as in your last example (defining the terms in the earlier ones is a good start).
> 1 votes

---
Tags: r, snp, gwas, variation
---


# Question
Title: What is the right way of calculating a Phred score by hand?
i am trying to calculate mean Phred scores for my sequencing data, but i feel not very comfortable about it. There are actually two ways of calculating. (I just use an existing sample) giving: 3 reads with a mean phred score of 20, 10 and 3 (meaning Error rates of 1%, 10% and 50%, respectively)

Way 1 to calculate: Just take the arithmetic mean of the q scores:

`(20+10+3)/3 = 11`

Way 2: Take the arithmetic mean of the error rates:

`(0,5+0,1+0,01)/3 ~= 0.2`

And then calculate the q score by

`-10*log10(0.2) = 6.99`

(https://gigabaseorgigabyte.wordpress.com/2017/06/26/averaging-basecall-quality-scores-the-right-way/)

The Author of this Blog claims, that Way 2 is the correct version, but i still do not get why.

Is somebody able to answer this question?

Thanks a lot

# Answer
Not sure if my explanation is any good but let's try...

First, let's convince ourselves that converting phred to probability is the right thing to do (*i.e.* your Way 2).

Imagine a read of length 10000 bp with half of the bases with Q=10 (or P=0.1) and the other half with Q=20 (or P=0.01). The number of wrong calls you expect is therefore 550:

```
n = 10000
(n * 0.01 * 0.5) + (n * 0.1 * 0.5) # -> 550 wrong calls

```

So the overall error rate is `550/10000 = 0.055` which corresponds to Q=12.6. If you follow your Way 2 you do get this answer:

```
((n * 0.5 * 0.01) + (n * 0.5 * 0.1))/n = 0.055
-10*log10(0.055) = 12.6

```

If you follow your Way 1 you get Q=15, which means an overoptimistic error rate of ~0.03 or ~316 wrong calls. Why is that?

Because phred is on a non-linear, logarithmic, scale: a unit change towards the high values is less important than a unit change towards the lower end. For example, across 10000 bases a change from Q=30 to Q=31 means a change from 10 to ~8 wrong calls - a difference of 2 only. But a change from Q=10 to Q=11 means 1000 to ~800 wrong calls: the difference is a substantial ~200 calls.

Having said this, the straight average of phred scores is not "wrong". It measures the overall phred score but it does not reflect the expected number of wrong calls, which is probably what you want in this case.
> 5 votes

---
Tags: ngs, reads, quality-control
---


# Question
Title: How to demultiplex a mix of single-indexed and dual-indexed samples
## The problem

If I have a sample sheet that contains both single-indexed and dual-indexed samples, I can split it up into two sample sheets and then run bcl2fastq on each one. However, when doing this, large Undetermined fastq files are generated. E.g., when processing the single-indexed samples, all the dual-indexed samples go to Undetermined. And when processing the dual-indexed samples, the single-indexed samples go to Undetermined.

Additionally, because they are being processed separately, if any single index A is part of a dual index A+B, then when processing the single indexes, an A+B may be mistaken for an A, so it would seem that they need to processed simultaneously to avoid this mis-assignment.

## The question

Given a sample sheet and directories of BCL files, how can such a set of sequencing data be demultiplexed correctly, either using bcl2fastq or the Picard tools?

To put it another way, I want to demultiplex a single sequencing run that contains both single-indexed and dual-indexed samples. It can be assumed that the indexes are sufficiently distinct such that any sample's index configuration is different from any other. But assuming that the different index configurations are not segregated to particular lanes of the sequencer, the question is how to demultiplex the files correctly such that both the single-indexed and dual-indexed samples are recognized.

## Attempts at a solution

### Using bcl2fastq directly

If a sample sheet of a mix of indexes is given to bcl2fastq (v2.17) directly, it produces the error

```
ERROR: bcl2fastq::common::Exception: Success (0): .../bcl2fastq2/src/cxx/lib/layout/BarcodeCollisionDetector.cpp(127):
Throw in function void bcl2fastq::layout::BarcodeCollisionDetector::validateNewBarcodeSizesAgainstExisting(const std::vector<long unsigned int>&) const
Dynamic exception type: boost::exception_detail::clone_impl<bcl2fastq::layout::BarcodeCollisionError>
std::exception::what: Barcodes have an unequal number of components.

```

> Barcodes have an unequal number of components.

### Using picard

It seems that it should be possible using the Picard tools, but I have not found a way to set up the inputs to ExtractIlluminaBarcodes and IlluminaBasecallsToFastq that processes this configuration correctly.

The syntax for supporting multiple index configurations is not entirely clear. But when using various combinations of N and '\*' on multiplex\_params.tsv and barcodes.txt required by the Picard tools, a large Undetermined file is still produced and actual sample fastq files are tiny, indicating it is not processing them correctly.

### Using multiple calls to bcl2fastq

As indicated in the discussion above, this suffers from the problem of having all the dual indexed samples going to the "Undetermined" file when processing the single-indexed samples, or vice versa, and creates two output directories of reports and stats, which must be merged.

### Padding empty indexes with N

By padding, I mean converting ACGTACGT to ACGTACGT+NNNNNNNN so that the single-index samples in the sample sheet are "dual" as well. This strategy seems that it would work except there is a bug in bcl2fastq that it treats "N" literally instead of as a wildcard. See the release notes for details.

# Answer
> 1 votes
I'm late to the party but simply run `cellranger mkfastq` twice, with different arguments for the mask and with `--filter-dual-index` for the double-indexed samples. In a second step you then have to disentangle the dual A+B reads from single-index A reads based on the readID's that are found in both of the demultiplexings. This is cumbersome and time+diskspace consuming, but it can be done. I would use some Unix tools for that but there are probably smarter ways.

EDIT: I mistakenly said `bcl2fastq`, this must be `cellranger mkfastq`. (thanks AmadeusDrZaius)

# Answer
> 0 votes
Have you just tried giving bcl2fastq one sample sheet with a mix of single and dual indices? I don't think what you are trying to do is a problem.

> if any single index A is part of a dual index A+B,

Well, it might be a problem if you did that. That was poor planning.

---
Tags: illumina, picard, demultiplexing, indexing
---


title: Is SBOL-VOWS operational?

text:
# Question
Title: Is SBOL-VOWS operational?
Synthetic Biology Open Language's Visual Ontology Webservice docs have example URLs but they all return a github pages 404.

EDIT: The SBOL-VOWS has site has been fixed but the section on the SBOL website still has links with a 404.

# Answer
> 3 votes
Many thanks for raising this issue. The web service is now working. Only the glyph related web service calls were affected due to restructuring folders at the SBOL-Visual GitHub repository. I now revived the glyphs using the sbol-visual-ontology GitHub repository. Over the summer, we will develop the web service further so that you can request different versions of glyphs. Currently, the web service works with the SBOL Visual v2.3 glyphs and we will provide support for SBOL Visual v3 glyphs soon.

The web service is available at: https://vows.sbolstandard.org/

Please let me know if there is anything I can clarify.

---
Tags: visualization, sbol, ontology
---


# Question
Title: Making Seurat object from a processed file
I have downloaded `log2(TPM/10+1)` values of `11,548 genes` and `9609 cells` from `GSE146026` (10x) in `tsv format` as the raw data is not available

I see patient IDs, cell barcodes, genes, even assigned clusters are here

```
> dput(head(a[1:5,1:5]))
structure(list(Cell_ID = c("10x_barcode", "patient", "time",
"sample_ID", "clst"), X10x_1 = c("10x_3288_t1_AAACATACCTTCCG-1",
"5", "1", "3288.1", "1"), X10x_2 = c("10x_3288_t1_AAACATACTCCTAT-1",
"5", "1", "3288.1", "1"), X10x_3 = c("10x_3288_t1_AAACATTGAACTGC-1",
"5", "1", "3288.1", "1"), X10x_4 = c("10x_3288_t1_AAACATTGCTGACA-1",
"5", "1", "3288.1", "2")), row.names = c(NA, 5L), class = "data.frame")
>

```

I want to make a Seurat object of that but really I don't know how to do that

I have tried this

```
    CreateSeuratObject(getwd(),gene.column = 2,cell.column = 1,unique.features = TRUE,strip.suffix = FALSE)

```

But get the following error:

```
Error in CreateAssayObject(counts = counts, min.cells = min.cells, min.features = min.features) :
  No cell names (colnames) names present in the input matrix

```

Anybody have ever dealt with such a case to help me please?

# Answer
The function you need is `CreateSeuratObject()` and not `Read10X()` as you start from TPMs.
> 2 votes

---
Tags: r, seurat, 10x-genomics
---


# Question
Title: How to introduce point mutations in a protein structure?
I want to introduce point mutations to a protein PDB file and automatically predict the 3D coordinates of a new amino acid residue in the mutated protein.

I used interactive Schrodinger's Maestro tool to mutate residues and minimize them to avoid clashes, however this time i want to automate the task and use some open source tools.

Is there any free software for this task?

# Answer
> 3 votes
**PDBFixer** and **pyMutate** are simple, lightweight, free and open-source. Both don't do any optimizing, though, so you may have to run some molecular dynamics afterward to get what would be the right "orientation".

**VMD** would work for sure, I believe **PyMOL** too. They are both graphical interfaces but you can start them in text only mode (or just `import pymol` in a python script and use the Python API). They are free for particulars and academics, but given their renown you may already be using them.

If you're looking for something a tad more complete and that may adapt the side-chain to the surrounding environment I'd go with **Modeller**, which provides two example scripts to do point mutations (here and there).

If you don't want to just get the side-chain of this amino-acid right but actually strive to see some general backbone adaptation too then use **PyRosetta**. It can probably do what all the others do, too, but given it is 1 GB big it may be overkill. Any of the package coupled with molecular dynamics would also do the trick.

---
Tags: protein-structure, point-mutation
---


# Question
Title: How to index residues by icode in biopython?
I am indexing a structure by the code below:

```
PDB_CODE = '1lve'
RES_CHAIN = 'A'
RES_NUM = 27

structure = PDBParser(is_pqr=True).get_structure(PDB_CODE, PDB_CODE+".pqr")
residue = structure[0][RES_CHAIN][RES_NUM]

```

However, the 1lve contains residues with the same RES\_NUM but different icodes 'A' to 'F' so my code either ignores them or throws an error if I pass '27A'. Does biopython allow this kind of structure indexing directly?

# Answer
Okay, I found the answer myself. It is possible to index not only by RES\_NUM, but by its id (check with residue.get\_id()), so the below example would work:

```
#structure[0][RES_CHAIN][het, resseq, icode]
structure[0][RES_CHAIN][' ', 27, 'A']

```
> 0 votes

---
Tags: python, protein-structure, pdb
---


# Question
Title: Adding legend to Botton and top annotation in complexheatmap
I have a heat map but top and bottom annotation don't have legend

```
Top_ha = HeatmapAnnotation(Group = anno_simple(top$Group,
                                               col =structure(hue_pal()(length(unique(top$Group))),
                                                              names = unique(top$Group)),
                                           border = TRUE,
                                           height = unit(5, "mm")),
                       "Total Mutation Burden"=anno_barplot(clin$TMB, height = unit(15, "mm"), width=unit(30, "mm"),gp = gpar(fill = structure(hue_pal()(length(unique(clin$TMB))),
                                                                                names = unique(clin$TMB)))))

Bottom_ha = HeatmapAnnotation(Response = anno_simple(clin$Response,
                                                 border = T,
                                                 col = structure(c("purple4", "grey"), names = c("Responders", "NonResponders")),
                                                 height = unit(5, "mm"),show_),
                          "TP53 Status" = anno_simple(clin$TP53_status,
                                                          border = TRUE,
                                                          col = structure(c("firebrick", "grey"), names = c("Mutant","Wild type")),
                                                          height = unit(5, "mm")))

Right_Anno <- HeatmapAnnotation(Signature = anno_text(right$Signature), which = "row")

Heatmap(mat,
        bottom_annotation = Bottom_ha,
        top_annotation = Top_ha,
        name = "Exposure",
        right_annotation = Right_Anno,
        show_row_dend = FALSE,
        show_parent_dend_line = FALSE,
        column_split = c(top$Group),
        column_title = c("Mutagenic", "DDR Impaired", "C>A/T Dominant"),
        col = colorRamp2(c( 0, 0.5), c("white", "red")),
        show_column_names = FALSE)

```

As you are seeing Response, group and total mutation burden don't have any legend and annotation

Please can you help me to put legend for annotation?

I have tried `annotation_lable` but says `Error: either provide` levels`or provide named`colors`.`

Also `show_annotation_name = TRUE,show_legend = TRUE` did not work

# Answer
*Answer from @ram-rs, converted from comment:*

For the `Group`, try using it without the `anno_simple`. Just `Top_ha = HeatmapAnnotation(Group=top$Group, "Total Mutation Burden" = ...)`, nothing else for `Group`. See if that generates a legend at all. We can customize later.hing else for Group. See if that generates a legend at all. We can customize later.
> 0 votes

---
Tags: r, heatmap
---


title: Interpreting this plot from GSEA

text:
# Question
Title: Interpreting this plot from GSEA
I have `RNA-seq` for two groups of patients: `Responders` to chemotherapy (n=9) versus `non-responders` to chemotherapy (n=24)

I have done `DESeq2` and obtain `fold change of responders versus non-responders` in which `positive fold change` means `down-regulation in responders` (if I am not wrong)

```
> head(df)
    gene       FC
1 286499 7.116071
2 260436 6.691189
3   2267 5.751322
4   5968 5.522079
5   4103 5.328008
6 442709 5.036087
>

```

I have done `GSEA` in `pre-ranked` way

got this plot for `KRAS pathway`

I interpret this as; `KRAS pathway` in being `down regulated` in `responders`

Am I right?

Please correct me

When I separated genes for up and down regulated

For up regulated genes I did not get any significant GSEA but for down regulated genes in responders I got this

SO does this mean that `KRAS pathway` is being `up regulated` in `responders` because I have used `down regulated genes in responders`?

# Answer
A *gene set enrichment analysis* (GSEA) tests for *enrichment* of a gene set within a ranked list of genes.

The primary outcome of the analysis is *enrichment* or *no enrichment*. Gene sets frequently include correlated genes, and that correlation can be positive (enhancers or co-expressors) or negative (e.g repressors). It's reasonable to include negatively-correlated genes within gene sets because some pathways involve both enhancement of a particular biochemical pathway, and suppression of a pathway that carries out opposite things.

I don't know about the KRAS pathway you're comparing against here, but a statement I would make based on interpreting that GSEA graph would be something like, "The KRAS.600\_UP.V1\_DN gene set is enriched when comparing the expression of responders vs non-responders."

I'm not quite sure how you separated up and down-regulated genes (code would be helpful with that). Your down regulated gene plot looks like it's ranked in the opposite direction to the combined plot, which is not what I expect. I think that determining an enrichment score separately for upregulated and downregulated genes (within your ranked list of genes) is a good idea, but it shouldn't change the direction of the enrichment curve within each category. It will, however, alter the form of the curve, by stretching or compressing values vertically and altering the start/end offset.
> 1 votes

# Answer
In Deseq2, it depends on contrast eg. if your `contrast <- c("Condition", firstC, SecondC)` is then `-ve` is downregulated in `FirstC` and same for `+ve`. It is always better to take look at count data for the confirmation (responders versus non-responders samples).

\[for more detail you can refer here\]\[1\]

\[1\]https://github.com/amarinderthind/RNA-seq-tutorial-for-gene-differential-expression-analysis/blob/master/RNA-Seq-DGE.R
> 1 votes

---
Tags: r, rna-seq, pathway, gsea
---


title: Selecting rows by partial match

text:
# Question
Title: Selecting rows by partial match
I have a data frame one column is sample name like

```
Sample

TCGA-VR-A8EX-01A-11D-A36I-01
TCGA-L7-A56G-01A-21D-A25X-01
TCGA-VR-A8ER-01A-11D-A36I-01
TCGA-JY-A6FB-01A-11D-A33D-01
TCGA-2H-A9GO-01A-11D-A37B-01
TCGA-Z6-AAPN-01A-11D-A402-01
TCGA-L5-A4OP-01A-11D-A25X-01
TCGA-L5-A8NV-01A-11D-A37B-01
TCGA-R6-A8W5-01B-11D-A37B-01
TCGA-S8-A6BW-01A-11D-A31T-01
TCGA-R6-A6Y0-01B-11D-A33D-01
TCGA-R6-A8WC-01A-11D-A37B-01
TCGA-S8-A6BV-01A-21D-A31T-01
TCGA-L5-A4OJ-01A-11D-A25X-01
TCGA-LN-A7HV-01A-21D-A350-01

```

I want to subset samples which have partial match with these name

```
TCGA-IC-A6RE
TCGA-IG-A4QS
TCGA-JY-A6F8
TCGA-JY-A6FB
TCGA-L5-A43E
TCGA-L5-A4OG
TCGA-L5-A4OH
TCGA-L5-A4OJ

```

I have used

```
data[grepl(c("TCGA-IC-A6RE", "TCGA-IG-A4QS", "TCGA-JY-A6F8", "TCGA-JY-A6FB", "TCGA-L5-A43E", "TCGA-L5-A4OG", "TCGA-L5-A4OH", "TCGA-L5-A4OJ", "TCGA-L5-A4ON", "TCGA-L5-A4OS", "TCGA-L5-A4OW", "TCGA-L5-A4OX", "TCGA-L5-A88T", "TCGA-Q9-A6FW"),data$Sample),]

```

But says that

```
Warning message:
In grepl(c("TCGA-IC-A6RE", "TCGA-IG-A4QS", "TCGA-JY-A6F8", "TCGA-JY-A6FB",  :
  argument 'pattern' has length > 1 and only the first element will be used

```

Likely only the first entry from the vector being grepped

Edit: As @StupidWolf suggested, I have tried:

```
> data[substr(data$Sample, 1, 12) %in%
    c("TCGA-IC-A6RE", "TCGA-IG-A4QS", "TCGA-JY-A6F8", "TCGA-JY-A6FB",
      "TCGA-L5-A43E", "TCGA-L5-A4OG", "TCGA-L5-A4OH", "TCGA-L5-A4OJ",
      "TCGA-L5-A4ON", "TCGA-L5-A4OS", "TCGA-L5-A4OW", "TCGA-L5-A4OX",
      "TCGA-L5-A88T", "TCGA-Q9-A6FW"),]

```

This works, thanks.

# Answer
> 0 votes
*Answer from @stupidwolf, converted from comment:*

You cannot use grepl with a vector, it must be a pattern.

I can suggest two things and you just put them together i hope, to know if something is found in a vector, you can use `data$Sample %in% c("TCGA-IC-A6RE", "TCGA-IG-A4QS")` etc ..

Of course `data$Sample` will not be in there, you need a bit of manipulation, so since your samples follow a certain pattern, you can use `substr(data$Sample,....)` . If you do `substr(data$Sample,....) %in% c(.....)` , this will give you the selection.

# Answer
> 2 votes
Not strictly bioinformatics, but a common problem:

`grep` assumes a regular expression, not a vector of patterns. Assuming your patterns are stored as a vector, so e.g.

```
patt <- c("pattern1", "pattern2", "pattern3")

```

you first have to write them as a regex:

```
patt.regex <- paste(patt, collapse = "|")

```

Output of `patt.regex` is then `"pattern1|pattern2|pattern3"` which is now ready for `grep` or `grepl`:

```
grep(patt.regex, data$Sample)

```

For exact matches it would be `paste(paste0("^", patt, "$"), collapse = "|")` to get `"^pattern1$|^pattern2$|^pattern3$"`.

# Answer
> 0 votes
When the string follows certain format there is no need for partial match, see this example:

```
# example data
data <- read.table(text = "
Sample
TCGA-VR-A8EX-01A-11D-A36I-01
TCGA-L7-A56G-01A-21D-A25X-01
TCGA-VR-A8ER-01A-11D-A36I-01
TCGA-JY-A6FB-01A-11D-A33D-01
", header = TRUE, stringsAsFactors = FALSE)

x <- c("TCGA-VR-A8EX", "TCGA-VR-A8ER")

data[ substr(data$Sample, 1, 12) %in% x, "Sample", drop = FALSE]
#                         Sample
# 1 TCGA-VR-A8EX-01A-11D-A36I-01
# 3 TCGA-VR-A8ER-01A-11D-A36I-01

```

Or when sample IDs lengths are not the same we can split on "-" then paste it back first 3 items, then subset as usual with `%in%`:

```
data[ sapply(strsplit(data$Sample, split = "-"),
             function(i) paste(i[1:3], collapse = "-")
             ) %in% x, "Sample", drop = FALSE]
#                         Sample
# 1 TCGA-VR-A8EX-01A-11D-A36I-01
# 3 TCGA-VR-A8ER-01A-11D-A36I-01

```

---
Tags: r
---


title: Permanent error in installing R package

text:
# Question
Title: Permanent error in installing R package
I want to install `NNLM`

```
install_github('linxihui/NNLM')

lang: error: unsupported option '-fopenmp'
make: *** [RcppExports.o] Error 1
ERROR: compilation failed for package ‘NNLM’
* removing ‘/Library/Frameworks/R.framework/Versions/3.6/Resources/library/NNLM’
Error: Failed to install 'NNLM' from GitHub:
  (converted from warning) installation of package ‘/var/folders/wk/qwgr9yqs6ts0czx__y_w8gch0000gq/T//RtmpA8PudK/filec5bce1f1309/NNLM_0.4.4.tar.gz’ had non-zero exit status
>
> install.packages('NNLM')
Warning in install.packages :
  package ‘NNLM’ is not available (for R version 3.6.3)
>

```

What should I do?

```
> session_info()
─ Session info ─────────────────────────────────────────────────────────────────
 setting  value
 version  R version 3.6.3 (2020-02-29)
 os       macOS Mojave 10.14.6
 system   x86_64, darwin15.6.0
 ui       RStudio
 language (EN)
 collate  en_GB.UTF-8
 ctype    en_GB.UTF-8
 tz       Europe/London
 date     2020-05-11

─ Packages ─────────────────────────────────────────────────────────────────────
 package      * version date       lib source
 abind          1.4-5   2016-07-21 [1] CRAN (R 3.6.0)
 assertthat     0.2.1   2019-03-21 [1] CRAN (R 3.6.0)
 backports      1.1.6   2020-04-05 [1] CRAN (R 3.6.3)
 broom          0.5.6   2020-04-20 [1] CRAN (R 3.6.3)
 callr          3.4.3   2020-03-28 [1] CRAN (R 3.6.2)
 car            3.0-7   2020-03-11 [1] CRAN (R 3.6.0)
 carData        3.0-3   2019-11-16 [1] CRAN (R 3.6.0)
 cellranger     1.1.0   2016-07-27 [1] CRAN (R 3.6.0)
 cli            2.0.2   2020-02-28 [1] CRAN (R 3.6.0)
 colorspace     1.4-1   2019-03-18 [1] CRAN (R 3.6.0)
 crayon         1.3.4   2017-09-16 [1] CRAN (R 3.6.0)
 curl           4.3     2019-12-02 [1] CRAN (R 3.6.0)
 data.table   * 1.12.8  2019-12-09 [1] CRAN (R 3.6.0)
 desc           1.2.0   2018-05-01 [1] CRAN (R 3.6.0)
 devtools     * 2.3.0   2020-04-10 [1] CRAN (R 3.6.3)
 digest         0.6.25  2020-02-23 [1] CRAN (R 3.6.0)
 dplyr          0.8.5   2020-03-07 [1] CRAN (R 3.6.0)
 ellipsis       0.3.0   2019-09-20 [1] CRAN (R 3.6.0)
 fansi          0.4.1   2020-01-08 [1] CRAN (R 3.6.0)
 farver         2.0.3   2020-01-16 [1] CRAN (R 3.6.0)
 forcats        0.5.0   2020-03-01 [1] CRAN (R 3.6.0)
 foreign        0.8-76  2020-03-03 [1] CRAN (R 3.6.0)
 fs             1.4.1   2020-04-04 [1] CRAN (R 3.6.3)
 generics       0.0.2   2018-11-29 [1] CRAN (R 3.6.0)
 ggplot2      * 3.3.0   2020-03-05 [1] CRAN (R 3.6.0)
 ggpubr       * 0.3.0   2020-05-04 [1] CRAN (R 3.6.3)
 ggsignif     * 0.6.0   2019-08-08 [1] CRAN (R 3.6.0)
 glue           1.4.0   2020-04-03 [1] CRAN (R 3.6.3)
 gtable         0.3.0   2019-03-25 [1] CRAN (R 3.6.0)
 haven          2.2.0   2019-11-08 [1] CRAN (R 3.6.0)
 hms            0.5.3   2020-01-08 [1] CRAN (R 3.6.0)
 labeling       0.3     2014-08-23 [1] CRAN (R 3.6.0)
 lattice        0.20-41 2020-04-02 [1] CRAN (R 3.6.3)
 lifecycle      0.2.0   2020-03-06 [1] CRAN (R 3.6.0)
 maftools     * 2.4.05  2020-05-10 [1] Github (PoisonAlien/maftools@6ddf526)
 magrittr       1.5     2014-11-22 [1] CRAN (R 3.6.0)
 Matrix         1.2-18  2019-11-27 [1] CRAN (R 3.6.3)
 memoise        1.1.0   2017-04-21 [1] CRAN (R 3.6.0)
 munsell        0.5.0   2018-06-12 [1] CRAN (R 3.6.0)
 nlme           3.1-147 2020-04-13 [1] CRAN (R 3.6.2)
 openxlsx       4.1.5   2020-05-06 [1] CRAN (R 3.6.2)
 pillar         1.4.4   2020-05-05 [1] CRAN (R 3.6.3)
 pkgbuild       1.0.8   2020-05-07 [1] CRAN (R 3.6.2)
 pkgconfig      2.0.3   2019-09-22 [1] CRAN (R 3.6.0)
 pkgload        1.0.2   2018-10-29 [1] CRAN (R 3.6.0)
 plyr           1.8.6   2020-03-03 [1] CRAN (R 3.6.0)
 prettyunits    1.1.1   2020-01-24 [1] CRAN (R 3.6.0)
 processx       3.4.2   2020-02-09 [1] CRAN (R 3.6.0)
 ps             1.3.3   2020-05-08 [1] CRAN (R 3.6.2)
 purrr          0.3.4   2020-04-17 [1] CRAN (R 3.6.2)
 R6             2.4.1   2019-11-12 [1] CRAN (R 3.6.0)
 RColorBrewer   1.1-2   2014-12-07 [1] CRAN (R 3.6.0)
 Rcpp           1.0.4.6 2020-04-09 [1] CRAN (R 3.6.3)
 readxl         1.3.1   2019-03-13 [1] CRAN (R 3.6.0)
 remotes      * 2.1.1   2020-02-15 [1] CRAN (R 3.6.0)
 reshape      * 0.8.8   2018-10-23 [1] CRAN (R 3.6.0)
 rio            0.5.16  2018-11-26 [1] CRAN (R 3.6.0)
 rlang          0.4.6   2020-05-02 [1] CRAN (R 3.6.2)
 rprojroot      1.3-2   2018-01-03 [1] CRAN (R 3.6.0)
 rstatix        0.5.0   2020-04-28 [1] CRAN (R 3.6.2)
 rstudioapi     0.11    2020-02-07 [1] CRAN (R 3.6.0)
 scales         1.1.0   2019-11-18 [1] CRAN (R 3.6.0)
 sessioninfo    1.1.1   2018-11-05 [1] CRAN (R 3.6.0)
 stringi        1.4.6   2020-02-17 [1] CRAN (R 3.6.0)
 survival       3.1-12  2020-04-10 [1] CRAN (R 3.6.3)
 testthat       2.3.2   2020-03-02 [1] CRAN (R 3.6.0)
 tibble         3.0.1   2020-04-20 [1] CRAN (R 3.6.3)
 tidyr          1.0.3   2020-05-07 [1] CRAN (R 3.6.2)
 tidyselect     1.0.0   2020-01-27 [1] CRAN (R 3.6.0)
 usethis      * 1.6.1   2020-04-29 [1] CRAN (R 3.6.2)
 vctrs          0.2.4   2020-03-10 [1] CRAN (R 3.6.0)
 withr          2.2.0   2020-04-20 [1] CRAN (R 3.6.2)
 zip            2.0.4   2019-09-01 [1] CRAN (R 3.6.0)

[1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library

```

Update: I tried installing llvm, but don't know what put here to finish installing llvm-openmp pc-133-235:build fi1d18

```
$ cmake ../llvm -DLLVM_ENABLE_PROJECTS=openmp -DCMAKE_C_COMPILER=<C compiler> -DCMAKE_CXX_COMPILER=<C++ compiler>
-bash: syntax error near unexpected token `newline'

```

# Answer
*Answer from @atpoint, converted from comment:*

There are dozens of threads in GitHub, SE and other support platforms that clearly explain what to do. (e.g. https://github.com/velocyto-team/velocyto.R/issues/2) The issue is that the compiler does not support openMP, so you have to install one that does.

*Answer from @devon-ryan, converted from comment:*

In conda on OSX we solve this by installing llvm-openmp, which provides the -fopenmp flag. Maybe the same is possible outside of conda (no clue). Don't compile it, find an easier way to install it.
> 0 votes

# Answer
try `devtools::install_github('linxihui/NNLM')`, if it doesn't work you need to install the gfortran binary on Mac.

If it still doesn't work, try to upgrade to 4.0.0.

When you have this kind of trouble, look at the issue page of the github page. You will almost always see that others have had the same issue.
> 1 votes

---
Tags: r, software-installation, software-usage
---


# Question
Title: Python package or CLI tool outputting mutations of sequence with respect to reference
I'm new to bioinformatics, and I'm looking for a tool that takes two FASTA files, one containing a reference, let's call it `A.fasta`, the other another sequence, let's call it `B.fasta`, and outputs just the mutations that `B` contains with respect to `A`.

I'm open as to the output format. Is there a standard for this?

Background: I'm trying to analyse Sars-CoV-2 sequences and identify mutations with respect to the reference. I've done a google search for a tool that does what I want to but couldn't find something straightaway. Is there no standard tool for this?

# Answer
If it's just a comparison of two sequences that you're looking for, then *diffseq* from the emboss toolkit should be okay:

https://www.bioinformatics.nl/cgi-bin/emboss/diffseq

I think this tool requires the sequences to be oriented in the same direction. If you discover a *lot* more variants than you expect, try reverse-complementing one of the sequences.
> 2 votes

---
Tags: sequence-alignment, fasta, samtools, covid-19, sars-cov-2
---


title: How does FDRtool work?

text:
# Question
Title: How does FDRtool work?
I have a question about using FDRtool. In the below code (on RNA seq data whose p values were acquired using Deseq2), the FDRtool was first used and thereafter p.adjust using the benjamini hochberg method. I didn't write this, its taken from someone else's code found online, in addition it is also the way the it is written in the pipeline at my university. Therefore, I assumed that this is how fdrtool is used and in my eyes it looks like they both adjust the p value with the same method. So that's why I asked, as I want to understand why you would need both these lines, which I've seen many people do:

```
FDR.ddsRes <- fdrtool(ddsRes$stat, statistic= "normal", plot = T)
ddsRes[,"padj"]  <- p.adjust(FDR.ddsRes$pval, method = "BH")

```

Dont they both correct for false discovery rate? Why do you need to run p.adjust after using fdrtool?

Thanks!!

# Answer
Something like this book chapter might help on understanding FDR. The packages or functions you call differ in their estimation of pi0 or the proportion of null hypothesis.

Basically, the `p.adjust(..method = "BH")` is a more conservative method with Benjamini-Hochberg, you are assuming the proportion of null features to be 1.

When you use fdrtool, you are estimating the proportion of null hypothesis and from there, the local FDR. This in theory gives you more power.

However, I have to note that in your code, you are not using `fdrtool` to estimate the fdr as explained above, because:

```
FDR.ddsRes <- fdrtool(ddsRes$stat, statistic= "normal", plot = FALSE)

```

With this, the raw p-values are stored in `FDR.ddsRes$pval` and the q value, which is the corrected p-value, is stored in `FDR.ddsRes$qval`. In your code you are basically taking the raw p-values, and applying Benjamin-Hochberg:

```
p.adjust(FDR.ddsRes$pval, method = "BH")

```

You should apply a two tailed test. From the wald stats, it should be:

```
derived_p = 2*pnorm(-abs(ddsRes$stat))

```

In your code, you are running it wrongly by passing the test statistic into `fdrtools`, we can compare the raw p-values:

```
head(data.frame(stat = ddsRes$stat,derived = derived_p,
DESeq2=ddsRes$pvalue,fdrtools=FDR.ddsRes$pval))

        stat    derived     DESeq2  fdrtools
1  1.2405796 0.21476110 0.21476110 0.1216358
2  1.2314421 0.21815755 0.21815755 0.1244054
3 -0.3292516 0.74196552 0.74196552 0.6811997
4 -0.4964574 0.61957174 0.61957174 0.5356140
5 -0.4101163 0.68172066 0.68172066 0.6088427
6 -1.8007943 0.07173531 0.07173531 0.0246428

```

To summarize, use the adjusted p-value returned by DESeq2. The code you have above is wrong unless there's a good reason to do a one-tail test.
> 1 votes

---
Tags: r, rna-seq, deseq2, rna, p-values
---


# Question
Title: What does presence/absence of black lines mean in IGV?
IGV (Integrated Genome Viewer) is a popular open-source tool for viewing alignment files. In my BAM file in IGV, some deletions have black lines running through them, and others don't. What causes this? Short black lines also appear adjacent to some mutations and insertions.

I generated the CIGAR for the SAM file myself, and am wondering if I made an error. I didn't find anything in IGV's documentation about this.

# Answer
Black line represents an indel. Read was lengthened to fit consensus.

No black line (blank), no read/sequence. Read likely ended or was truncated.
> 1 votes

---
Tags: sam, igv, cigar
---


# Question
Title: Retain only part of a file name and fasta header in fasta directory
I have a `MainDirectory/hundreds_of_subdirectories/thousands_of_fastas` structure. For each fasta in the bottom level, I am trying to change the file name as well as the first line in the file to retain only some parts of the file.

So that the following text, which makes up both the file name and the first name of the file:

```
947-length-1150-cov-1000|contig:JAECWU010000213.1|slice:817050-818200|uce:uce-1452|match:817550-817700|orient:{'+'}|probes:4.unaligned.fasta

```

becomes just:

```
uce-1452.unaligned.fasta

```

For editing the first line in the file, I have tried

```
sed '1 s/^[^|uce-]*(|uce-[0-9]).*/\1/' hundreds_of_subdirectories/*

```

but this isn't the ticket..

Then messing around with some kind of loop...

```
for i in *\*.unaligned.fasta
do
  sed -E 's/^[^|uce-]*(|uce-[0-9]).*/\1/'
done

```

For the file name, it was suggested that I use rename instead of sed, but with sed I have tried something like:

```
for x in hundreds_of_subdirectories/thousands_of_fastas*unaligned.fasta
do
  echo $x | sed -r 's/^[^|uce-]*(|uce-[0-9]).*/mv' -v "\0" "|uce-[0-9].*/\1/.unaligned.fasta"/
done

```

Am I even barking up the right tree? Thanks!

# Answer
I think you're working too hard to capture the `uce:` expression.

Try

```
sed -Ee '1s/^.+uce:(uce-[0-9]+).+/\1.unaligned.fasta/' file_glob

```

to replace the first line content. You may want to add the in-place flag after testing the expression on a few files.

For the `mv`, use command substitution, which runs the command in a sub-shell and uses the output value as the parameter:

```
for f in file_glob
do
    mv -- "$f" "$(echo \"$f\" | sed '...')"
done

```
> 1 votes

# Answer
If you are using `bash`, you can enable the `globstar` option (`shopt -s globstar`) which allows `**` to match recursively. That way, `**/*fasta` will find *all* files whose name ends with `fastq` in this directory and all of its subdirectories. So that will let you easily iterate over all the target names. Next, if you have the perl `rename` command (called `rename` on Debian and Ubuntu and similar systems, `prename` or `perl-rename` on others), you can do it quite easily.

Note that `perl-rename` has the `-n` option that makes it print out the renaming it would do without actually renaming anything. So I suggest you first run this to test what you will be doing:

```
shopt -s globstar
for file in **/*.fasta; do
   echo perl -i.bak -pe 's/^\s*>.*uce:(uce-\d+).*(\.unaligned.fasta)/$1$2/' "$file"
   rename -n 's/^.*uce:(uce-\d+).*(\.unaligned.fasta)/$1$2/' "$file"
done

```

That will `echo` the perl command, so you can try it manually first and make sure it works, and will print out the `rename` operations without actually doing anything. If you are satisfied it will work, run again without the `echo` or `-n`:

```
for file in **/*fasta; do
    perl -i.bak -pe 's/^\s*>.*uce:(uce-\d+).*(\.unaligned.fasta)/$1$2/' "$file
    rename 's/^.*uce:(uce-\d+).*(\.unaligned.fasta)/$1$2/' "$file"
done

```

Or, if you don't have `perl-rename`, you can use the shell's string manipulation for this:

```
for file in **/*.fasta; do
    perl -i.bak -pe 's/^\s*>.*uce:(uce-\d+).*(\.unaligned.fasta)/$1$2/' "$file
    tmpFileName="${file#*uce:}"
    mv -- "$file" "${tmpFileName/|*/.unaligned.fasta}"
done

```

### Explanation

* `shopt -s globstar`: this enables the `globstar` bash option for the current shell session. It only needs to be run once, its effects will last until you close the terminal or start a new shell session. As for what it does, see `man bash`:

> globstar
>
> If set, the pattern `**` used in a pathname expansion context will match all files and zero or more directories and subdirectories. If the pattern is followed by a `/`, only directories and subdirectories match.
* `for file in **/*fasta; do ... ; done`: iterate over all files (and directories) whose name ends in `.fasta`, saving each file name as `$file`.
* `perl -i.bak -pe 's/^\s*>.*(uce-\d+).*(\.unaligned.fasta)/$1$2/' "$file`: this changes the fasta header. The `-i.bak` option tells `perl` to edit the file in place, but to save the original file with a `.bak` extension. This can help you revert any changes if there is a mistake. Once you have finished, you could run `rm **/*fasta.bak` to delete the backups. Alternatively, just remove the `-i.bak` if you don't want to keep backups.

The `-p` option tells `perl` to process its input line by line and print every line after applying whatever script you give it with `-e`.

The `s/old/new/` is the substitution operator, same as `sed`. The regular expression is looking for 0 or more whitespace characters and then a `>` at the beginning of the line (`^\s*>`). It then allows for 0 or more characters (`.*`) until it finds the string `uce:` followed by `uce-` and one or more digits (`uce:(uce-\d+)`). The parentheses "capture" the pattern, storing it as `$1`. Next, we want 0 or more characters again until we reach `.unaligned.fasta` (`.*(\.unaligned.fasta)`). We capture the `.unaligned.fasta` as `$2` (second set of parentheses) only to avoid typing it out again. Note that this assumes you will never have anything after the `.unaligned.fasta`. Finally, we replace everything matched above with just `$1` (the `uce-numbers`) and `$2` (`.unaligned.fasta`).
* `rename 's/^.*uce:(uce-\d+).*(\.unaligned.fasta)/$1$2/' "$file"`: the rename opration. This will match everything from the beginning of the file name until the last occurrence of `uce:uce-NNNN`, capture the `uce-NNNN` as `$1` then skip the remaining characters until `.unaligned.fasta` and finally rename as just `$1$2` (the capturing logic is the same as above).
* `tmpFileName="${file#*uce:}"`: The `${var#pattern}` syntax means "print the value of `$var`, but remove the shortest match of `pattern` from the front of the variable's value". The `${var//pattern/replacement}` syntax means "print the value of `$var` and replace the first occurrence of `pattern` with `replacement`. Here, the pattern is `|*` which means "anything after a `|`" (note that this is shell stuff and is using globs and not regular expressions). This is probably easier to understand with an example:

```
$ echo "$file"
947-length-1150-cov-1000|contig:JAECWU010000213.1|slice:817050-818200|uce:uce-1452|match:817550-817700|orient:{'+'}|probes:4.unaligned.fasta
$ tmpFileName="${file#*uce:}"
$ echo "$tmpFileName"
uce-1452|match:817550-817700|orient:{'+'}|probes:4.unaligned.fasta
$ echo "${tmpFileName/|*/.unaligned.fasta}"
uce-1452.unaligned.fasta

```

So, given `tmpFileName="${file#*uce:}"`, `"${tmpFileName/|*/.unaligned.fasta}"` is the desired file name and this is what we pass to `mv`.
> 1 votes

---
Tags: fasta, bash
---


# Question
Title: RDAVIDWebService (R package) manual installation, or alternative tool
Recently, I had to re-install all my R packages (R version 4.1.0 (2021-05-18) -- "Camp Pontanezen") and it turns out that `"RDAVIDWebService"` is not available anymore on Bioconductor:

> `package ‘RDAVIDWebService’ is not available for Bioconductor version '3.13'`

Is there an another way to install it (couldn't find the files on the Bioconductor site to install the package manually) or is there another R package that allows to use the DAVID tool?

# Answer
> 1 votes
That package was deprecated in Bioconductor version 3.13 after several attempts to contact the package author/maintainer.

Your options include reverting back to Bioconductor version 3.12 or using one of the other options for GO/pathway enrichment analyses, e.g. topGO, enrichR, or clusterProfiler.

---
Tags: r, bioconductor, gene-ontology
---


title: Solutions to reference bias issue

text:
# Question
Title: Solutions to reference bias issue
I've recently learned about the reference bias issue - inability to properly map NGS reads in certain genomic regions caused by the fact that our standard genomic references are linear and do not fully represent the genomic variation present in the population. According to what I read this has multiple consequences: missing structural variants or variants in highly variable regions like HLA for humans.

Now I wonder what is the best way to address this. Of course, "best" is context-dependent but I don't have any specific context so general pros and cons will be highly appreciated. I know about two approaches but probably there are more. The first one is an alt-aware mapping that's used in DRAGEN by Illumina. As far as I understand, they just keep multiple linear alternatives for some regions and try doing a "smarter" mapping for those regions. The second approach that I've found is using graph references where the reference is not linear anymore but rather a graph and mapping is taking into account this graph structure.

Thanks in advance.

# Answer
*Answer from @throckmorton, converted from comment:*

Here is a recent approach to this problem, using multiple linear references:

https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-02229-3
> 2 votes

# Answer
Slightly older than @gringer 's reference, but still from Richard Durbin himself

"Variation graph toolkit improves read mapping by representing genetic variation in the reference"

https://pubmed.ncbi.nlm.nih.gov/30125266/
> 2 votes

---
Tags: variant-calling, read-mapping, reference-genome, structural-variation
---


# Question
Title: RNA_Seq data aligned used uniquely or multi mapped reads impact on result interpretation
I have some transcriptomic (Whole) sequencing data that I should analyse. I would like to do raw data alignment to a reference genome taking into account the multi mapped reads and uniquely mapped reads. I would like to know the consequences on my transcriptomic analysis results : 1. if I should only focus on uniqueMapping reads (mapped to one loci) and neglect MultiMapping read (mapped to several loci or variants)! 2. or take them all together

So my questions are :

1- Considering unique and multi-aligned reads, wouldn't this be a source of bias in my results interpretation?

2- By hiding the multi-mapped reads, could I lose information on the expression of important genes?

thank you in advance

# Answer
> 3 votes
There are several methods/tools that take into account multi-mappers and deal with them in reasonable ways (e.g. pseudo-alignment methods like salmon and more traditional alignment-based methods like RSEM). If you completely ignore multi-mappers, then yes, you will be losing information, which may or may not be valuable to you.

Both of those methods utilize the expectation-maximization algorithm (or a variant thereof) to determine the relative contributions of multi-mapping reads to the various mapping locations. With salmon, you can forego alignment completely and end up with very accurate gene quantifications in hand in like 20 minutes if you want.

---
Tags: rna-seq, sam, ngs, genomics, transcriptome
---


# Question
Title: Bruker MALDI-TOF bacteria species identification scoring algorithm
Just wondering whether anyone can point me to a research paper which describes how the scoring values are generated by the Bruker software for species identification of bacteria in MALDI-TOF MS. For example, there are countless papers describing the Bruker output score as green (species-level confidence), yellow (genus-level), or red (neither species- nor genus-level confidence). However, none of the papers that I see are discussing how these numbers are actually calculated within the Bruker system. Even the Bruker brochure glosses over the point: "After the acquisition of the spectral data has been completed, a report is generated. The result for each sample is clearly listed under ‘Organism (best match)’ accompanied by the resulting score and an appropriate ‘traffic light’ color scheme." Can anyone point me in the direction of a paper or at least some online information such as wikipedia describing how Bruker came up with these numbers? It's not just an arbitrary number, right?

Another generic description of the scores is in this presentation: "Higher the log (score), higher the similarity between mass spectrum of isolate & the database entry in the reference library."

# Answer
> 1 votes
> For identification scoring, the software correlates signal intensities of matched signals of mass spectra. The three scores obtained from such a procedure are multiplied and normalized to a value of 1,000 and its common logarithm. Log scores over 2 are considered as reliable identification of a bacterial species (Table 4), whereas log scores over 1.7 generally indicate reliable identification of bacterial genera. Log scores of 3 are obtained when spectra are matched with themselves.

Sauer, S., Freiwald, A., Maier, T., Kube, M., Reinhardt, R., Kostrzewa, M., & Geider, K. (2008). Classification and identification of bacteria by mass spectrometry and computational analysis. *PloS one, 3*(7), e2843.

Going further, from a random thread on ResearchGate:

https://www.researchgate.net/post/Does\_anyone\_know\_how\_is\_the\_p-value\_calculated\_by\_Bruker\_Biotyper\_during\_MALDI-MS\_fingerprinting\_of\_bacteria

> Number of peaks matched in your unknown sample divided by the number of peaks in the unknown sample, this is given a score (out of 10)
>
> Number of peaks matched in your unknown sample divided by the number of peaks in the reference sample, this is given a score (out of 10)
>
> Score of relative intensities of matching peaks (out of 10)
>
> Total score out of 1000
>
> Log10 score out of 3
>
> There is no relationship of the score with a p-value nor you can calculate the statistics because the bruker database is not open and you can not know the values for what the score was calculated to do the statistical analysis

---
Tags: bacteria, mass-spectrometry
---


# Question
Title: Construct the Overlap Graph of a Collection of k-mers
I am doing a course on graph algorithms in genome sequencing and the current assignment involves building an overlap graph from a bunch of k-mers. Formally the problem is as follows:
**Input:** A collection of Patterns of k-mers.
**Output:** The overlap graph in the form of an adjacency list.
I have a working algorithm that handles most cases except for self-loops. The problem is that I don't know when/why I should be adding self-loops. Any help on this would be greatly appreciated. Below are some examples illustrating my question and my current solution. Thank you!

*Example 1*:
Input:

ACT
CTT
TTT

Output:

TTT-\>TTT
ACT-\>CTT
CTT-\>TTT

*Example 2*:
Input:

CCCC
Output:

CCCC-\>CCCC

*Example 3*:
Input:

CT
TG
TG
TC
TT
TC
Output:

TC-\>CT
CT-\>TC,TG,TT
TT-\>TC,TG,TT

The goal here is to represent a collection of kmers as an overlap graph where unique kmer will be represented as a node in the graph. Directed edges will connect nodes whose last $k-1$ characters overlap with the other node's first $k-1$ characters. Therefore, repeated kmers that overlap with a other kmer should be reflected as self loops. That's what I understood, however I am looking for a more general rule/solution.

Here is another example to illustrate this.

Input:

```
CT
TT
TT
TT
TT
TT

```

Output:

```
CT->TT
TT->TT

```

Current code:

```
/*
 * Creates an overlap graph from a list of k-mers
 *
 * \param patterns      a vector of string k-mers
 *
 * \return a string containing an adjacency list representation of the
 * overlap graph as described in the problem specification
 *
 * Output Format. The overlap graph Overlap(Patterns), in the form of an
 * adjacency list. Specifically, each line of the output represents all
 * edges (u,v) leaving a given node u in the format “u -> v”, where u and v
 * are both k-mers. If a given node u has multiple edges leaving it (e.g. v and w),
 * the destination nodes are comma-separated in any order. For example, if
 * there exist nodes ACG, CGT, and CGG, the resulting line of the adjacency
 * list would be: ACG -> CGT,CGG
 */

#include <algorithm>
#include <iostream>
#include <map>
#include <vector>
#include <set>
#include <string>

std::string overlap_graph(const std::vector<std::string>& patterns) {
    // Overlap Graph of the k-mers collection patterns, in which a vertex exists for
    // each unique k-mer, and an edge (u,v) exists if the last k-1 characters of u are
    // equal to the first k-1 characters of v.

    if (patterns.size() == 1) return patterns[0] + "->" + patterns[0];
    // Build adjancency list O(n^2)
    std::map <std::string, std::vector<std::string>> adj;
    for (const auto& kmer_u : patterns) {
        // Add kmer as key if it doesn't exist in the map and find its neighbours.
        if (adj.find(kmer_u) == adj.end()) {
            adj[kmer_u] = {};
            for (const auto& kmer_v : patterns) {
                if ((kmer_u != kmer_v) &&
                    (kmer_u.substr(1, kmer_u.length() - 1) == kmer_v.substr(0, kmer_v.length() - 1)) &&
                    (std::find(adj[kmer_u].begin(), adj[kmer_u].end(), kmer_v) == adj[kmer_u].end())) {
                    adj[kmer_u].push_back(kmer_v);
                }
            }
        } // Otherwise continue.
    }

    // Convert to string O(n).
    std::string ans;
    for (auto const& key_value : adj) {
        std::string kmer_u{ key_value.first };
        std::vector<std::string> kmers_v{ key_value.second };
        if (!kmers_v.empty()) {
            ans.append(kmer_u + "->" + kmers_v[0]);
            for (std::size_t kmer_v{ 1 }; kmer_v < kmers_v.size(); kmer_v++) {
                ans.append(',' + kmers_v[kmer_v]);
            }
            ans.push_back('\n');
        }
    }
    return ans;
}

void tests() {
    /* Sample Dataset */
    std::vector<std::string> patterns{ {"AAG"}, {"AGA"}, {"ATT"}, {"CTA"}, {"CTC"}, {"GAT"},
                                       {"TAC"}, {"TCT"}, {"TCT"}, {"TTC"} };
    std::string answer{ "TTC->TCT\nCTA->TAC\nAAG->AGA\nCTC->TCT\nTCT->CTA,CTC\nAGA->GAT\nGAT->ATT\nATT->TTC\n" };
    std::cout << answer << std::endl;
    // std::cout << (overlap_graph(patterns) == answer) << std::endl;
    std::cout << overlap_graph(patterns) << std::endl;
}

void solution() {
    std::vector<std::string> patterns;

    std::string pattern;
    while (std::cin >> pattern) {
        patterns.push_back(pattern);
    }
    std::cout << overlap_graph(patterns) << std::endl;
}

int main() {
    //tests();
    solution();
}

```

# Answer
You will want to create a self loop for a $k$-mer $s$ whenever $s\_1,...,s\_{k-1} = s\_2,...,s\_k$. This will only happen when all of the characters of the $k$-mer are identical, as $s\_1=s\_2, s\_2=s\_3,...,s\_{k-1}=s\_k$.
> 1 votes

# Answer
The usual way to avoid loops in linked lists is to keep a record of the things that have already been used / visited, and stop when an already visited thing has been seen again (or when there is nothing more to visit).
> 0 votes

---
Tags: ngs, k-mer, homework, c++, graphs
---


# Question
Title: Smith-Waterman traceback with large gaps
I am currently trying to understand how the traceback algorithm is supposed to work for the smith-waterman algorithm as my current understanding breaks down in case of a large alignment gap.

Assume the sequences

```
3GA12CA
GA56CA

```

The obvious best alignment would be (or some other variation thereof)

```
3GA12-CA
-GA--6CA

```

However, assuming an affine gap penalty function of `x => x + 2` and a binary similarity score (i.e., if equal then 1 else -1) would lead to a scoring matrix

```
    3 G A 1 2 C A
  0 0 0 0 0 0 0 0
G 0 0 1 0 0 0 0 1
A 0 0 0 2 0 0 0 0
6 0 0 0 0 0 0 0 0
C 0 0 0 0 0 0 1 0
A 0 0 0 1 0 0 0 2

```

When doing the traceback I would start with the largest value towards the lower right corner of the matrix

1. (6,6)
2. (5,5)
3. Now I do not know how to find the next best alignment (in this example at (3,2)) as neither going up nor going to the left nor going diagonal is going to lead to a score greater 0 meaning I have no next best value to go on.

How is the traceback function supposed to continue here?

Additionally, I read that the traceback should start at the highest score in the matrix itself. However, if we would extend the first subsequence match length by one, this would mean we miss the right alignment option completely.

# Answer
This question is based on a misunderstanding. The Smith-Waterman algorithm calculates the best `local` segment. You have to start at the $max(H\_{i,j})$ point and run the traceback as described. It will (most likely) return only a local segment that, however, has the maximum value in terms of alignment score.

To get a `global` alignment, another algorithm has to be used. Needleman-Wunsch would be an example.

The difference between `local` and `global` is that a `local` alignment only alignes a subsequence of the provided sequences whereas a `global` alignment gives the best alignment of the whole sequences.
> 1 votes

---
Tags: sequence-alignment
---


# Question
Title: Access base aligned to particular reference position
**The short version**: If I have a SAM record, is there any simple way to retrieve the base aligned to a particular reference position without computing a pileup?

**The long version**: I'm using pysam to write some genotyping code. I have a BAM file with reads aligned to some amplicons. I'd now like to genotype a couple of positions of interest in each read while maintaining the phase information. For example, at amplicon 1, if read 1 has `A` for SNP1 and `T` at SNP2, then I want to report an `A,T` genotype for read 1 at amplicon 1.

The typical way to genotype SNPs of interest is to compute a pileup and look at all of the aligned bases position-by-position. The pysam API includes the following example of how to work with pileups.

```
import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb" )
for pileupcolumn in samfile.pileup("chr1", 100, 120):
    print ("\ncoverage at base %s = %s" %
           (pileupcolumn.pos, pileupcolumn.n))
    for pileupread in pileupcolumn.pileups:
        if not pileupread.is_del and not pileupread.is_refskip:
            # query position is None if is_del or is_refskip is set.
            print ('\tbase in read %s = %s' %
                  (pileupread.alignment.query_name,
                   pileupread.alignment.query_sequence[pileupread.query_position]))

samfile.close()

```

However, I lose the phasing information if I analyze the data position-by-position. Is there any way to flip this example so that I'm iterating over reads in the outer loop, and then over positions in the inner loop? In other words, **if I'm iterating over reads, is there any simple way to retrieve the base aligned to a particular reference position *without* computing a pileup**?

As a potential solution, I've though about storing genotype information in a nested dictionary indexed first by read ID and then by position, but depending on coverage that could be incredibly memory intensive.

# Answer
> 4 votes
You can use the `get_reference_positions()` function in pysam to get a vector of the positions. Searching that for your positions will allow you to output phased calls.

```
for read in bam.fetch(some_position):
    positions = read.get_reference_positions():
    if some_other_position in positions:
        idx = positions.index(some_position)
        idx2 = positions.index(some_other_positions)
        ... and so on...

```

That's not exactly correct syntax, but you get the idea.

# Answer
> 3 votes
I had the same problem as you, and so I wrote my own code and sharing it here. No guarantee that it's bug-free ;-)

This code takes advantage of the cigar string. It returns None if the base has been deleted or clipped, or if the bam file did not store the sequence.

```
from typing import Optional
import pysam

def find_base_in_alignment(alignment: pysam.AlignedSegment,
                           pos: int,
                           bam_stores_revcomp: bool = False) \
        -> Optional[str]:
    idx_q = 0
    idx_r = pos - alignment.reference_start
    if bam_stores_revcomp:
        seq = alignment.query_sequence
    else:
        seq = alignment.get_forward_sequence()

    if seq is None:
        return None

    for op, l in alignment.cigartuples:
        ref_consumed = op in {0, 2, 3, 7, 8}
        query_consumed = op in {0, 1, 4, 7, 8}

        if ref_consumed:
            idx_r -= l
        if query_consumed:
            idx_q += l

        if idx_r < 0:
            if query_consumed:
                # base is in query between idx_q-l , idx_q
                base = seq[idx_q + idx_r - 1]
                return base
            else:
                # position has been deleted
                return None

```

Let's say you want to genotype all reads by looking at the base at chromosome `chrom` and position `pos`:

```
bamfile = pysam.AlignmentFile(filename, "r" )
for alignment in bamfile.fetch(chrom, pos-1, pos+1):
    base = find_base_in_alignment(alignment, pos)
    # ... do something with it

```

Note that the parameter `bam_stores_revcomp` is important. Some aligners will store the original query sequence in the bam, and other aligners (like minimap2) will store the forward sequence (i.e. the reverse complement whenever a read was aligned to the reverse strand). If you use the former, set `bam_stores_revcomp` to `False`, otherwise set it to `True`

# Answer
> 1 votes
This is not \[yet\] a complete answer, but hopefully it helps get you on the right path.

I wrote some code to retrieve the subsequence from reads in a BAM file that are actually mapped to the reference, which required calculating mapped base locations. Here's a subsection of the code that does this (slightly paraphrased):

```
  my $refPos = $F[3];
  my $cigar = $F[5];
  $cigar =~ s/[0-9]S$//;
  my $seq = $F[9];
  my $qual = $F[10];
  my $startTrim = 0;
  my $matchLen = 0;
  while($cigar =~ s/^([0-9]+)([MIDNSHP=X])//){
    my $subLen = $1;
    my $op = $2;
    if($op eq "S"){
        $seq = substr($seq, $subLen);
        $qual = substr($qual, $subLen);
      }
    }
    if($op =~ /[M=XI]/){
      $matchLen += $subLen;
    }
  }
  $seq = substr($seq, 0, $matchLen + $addSeq);
  $qual = substr($qual, 0, $matchLen + $addSeq);

```

So, some things to take note of:

* The entire CIGAR string (excluding terminal soft clips) should be parsed up to the required base location
* Any soft clips at the start of the sequence need to be accounted for
* Only update the read location counter if there's a match, mismatch, or insertion

---
Tags: sam, pysam
---


# Question
Title: Detecting SARS-CoV-2 variants from the mixed virus population
I have a fastq file from Nanopore sequencing data that contains reads from **both** the UK and South Africa variants of SARS-CoV-2.

The variants are identified by three key mutations in the receptor-binding domain of Spike.

**Aim** Firstly, I want to determine the number of reads then belongs to UK and South Africa variant. Then, I will create a bam file from the fastq file using the SARS-CoV-2 reference genome.

**Questions**

1. How do I determine the reads belonging to each variants?
2. Are there any tools available to classify the reads based on the pre-defined mutation criteria?

# Answer
> 2 votes
Based on a quick request from someone else, I wrote some code a couple of years ago for assigning read groups based on variants at particular locations; this was set up for HIV variant splitting, so I expect it should also work on SARS-CoV-2 nanopore reads as well:

https://gitlab.com/gringer/bioinfscripts/-/blob/master/samVarSplitter.pl

It processes a SAM stream \[i.e. in text mode\] using Perl, so might be a bit slow, but should be fine for a few hundred thousand sequences. You'd use it as follows:

```
samtools view -h mapped_reads.bam | \
  samVarSplitter.pl [-ref <refName>] [-pos <int>] | \
  samtools sort > grouped_reads.bam

```

Then the reads can be split into different files based on their group using `samtools split`:

```
Usage: samtools split [-u <unaccounted.bam>] [-h <unaccounted_header.sam>]
                      [-f <format_string>] [-v] <merged.bam>

```

---
Tags: sam, variant-calling, nanopore, covid-19, sars-cov-2
---


# Question
Title: Annotating a .vcf with centimorgan information
Some programs (e.g. shapeit4) automatically annotate an `INFO` tag into a .vcf file which gives the cumulative genetic distance in cM between each SNP:

```
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=10/11/2020 - 21:29:39
##source=G2H
##contig=<ID=9>
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=AC,Number=1,Type=Integer,Description="Allele count">
##INFO=<ID=CM,Number=A,Type=Float,Description="Interpolated cM position">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Phased genotypes">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##bcftools_viewVersion=1.9-207-g2299ab6+htslib-1.9-271-g6738132
##bcftools_viewCommand=view -s -1_-1 ukb_cal_chr9_v2.HumanOrigins.common_snps.HumanOriginsInds.fiftyPercentEthnicUKBInds.PHASED.bgz.vcf.gz; Date=Mon May 24 14:15:51 2021
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  -1_-1
9   257034  rs573860    A   G   .   .   AC=0;AF=0.388699;CM=0;AN=2  GT  0|0
9   270321  rs17720310  A   G   .   .   AC=0;AF=0.0820595;CM=0.007994;AN=2  GT  0|0

```

However, this program is really used for phasing not for .vcf annotation.

Ideally I'd like a simple command line program which takes input the .vcf and a genetic map and interpolates the genetic distance and adds a cumulative map. This could be done manually with a custom script, but quite frankly I can't be bothered to do that.

Does anyone know of a tool that exists which will do this?

I posted this to biostars a while ago and received one reply which linked to this post, which includes links to tables containing cM locations for different rs numbers.

EDIT: I am interested in doing this for human data (although the same would be the case for any diploid organism).

Sample vcf:

https://pastebin.com/raw/YSbuU0qM

# Answer
> 2 votes
I think a useful tool to address your question is QCTOOL v2. Particularly, the option you are looking for is -annotate-genetic-map. This tool has been developed originally for GWAS datasets. It has multiple inputs data format.

If you could provide input data as an example, I will try working on it to provide the necessary steps to be done.

Hope this answers your question.

---
Tags: vcf, sequence-annotation
---


# Question
Title: Remove repetitive region vcf file using repeatmasker bed file
I have a 1000 genomes vcf file for chromosome 14. And I want to remove variants in repetitive regions flagged by repeatMasker. I have a bed file for those repeats. I downloaded that file from UCSC. How can I eliminate/remove those sites in vcf file using repeat masker's bed file? Insights will be appreciated.

# Answer
> 2 votes
A couple fast alternatives:

```
$ bedops --chrom chr14 -n 1 <(vcf2bed < variants.vcf) <(sort-bed repeats.bed) > answer.bed
$ bedops --chrom chr14 -n 1 <(vcf2bed < variants.vcf) <(rmsk2bed < repeats.rmsk) > answer.bed

```

Using `--chrom chr14` focuses set operations on that chromosome.

# Answer
> 0 votes
`vcftools` with the `--exclude-bed` option or `bedtools` with the `subtract` command should do the trick.

---
Tags: vcf, snp, bedtools, repeatmasker, variants
---


title: where to find sample contig data

text:
# Question
Title: where to find sample contig data
Where or How can I find contigs? I am trying to learn Bioinformatics and I want to make a reference-based gene search. I also want to align contigs with a given reference sequence. From NCBI other databases I can find only reads and complete genome sequences. Please tell me how can I find contig sequences of a species?

# Answer
> 2 votes
Most of the reference genomes on NCBI are made of contigs or scaffolds (which is just a fancy name for contigs glued together by a bunch of unknown nucleotides). Hence, by download the whole genome, you also download all the contigs.

If you search for a specific contig (presumably a known locus) of a specific species, you can search NCBI for the species name + locus / scaffold / contig name and if it was sequenced it will pop up in the "genes" section of NCBI. For example, here is a search for *Pan troglodites* contig NC\_036898.1 (I typed literally Pan troglodites NC\_036898.1 in the search bar). You can even just type a species name, and click on "genes" and see loads of all the individual contigs sequenced.

--- edit ---

There is a bunch of ways to generate something we could call contig. However, usually what we mean is "continuous sequence assembled from sequencing reads" and it's generated from whole-genome sequencing reads. We sequence a genome, assemble the reads and what we get is a set of contigs as a representation of a haploid reference genome. Some of them will be short fragments of nucleotides with unknown functions, some of them can be long and contain hundreds of genes. The important part is that it's a continuous sequence.

With a bit more relaxed definition of "contig", we could also accept as contigs outputs of targeted sequencing (including old-fashioned sanger). It's a continuous stretch of nucleotides, usually bordered by a pair of primers. I am not 100% sure, but my best guess that most of the "gene" sequences in NCBI will be sanger-based (but I should disclaim I never investigated too much where all the NCBI sequences come from).

---
Tags: sequence-alignment, gene, assembly
---


# Question
Title: Important genes beyond PAM50 for breast cancer classification
*Note: this question has also been asked on Biostars*

I am currently trying to complete a Breast Cancer Classification task using Neural Networks. I have experimented with using my full dataset of gene expressions(1000+) and managed to improve my results by using only the PAM50 genes. I was wondering though what other genes I could add to the list to even better improve my results.

I have tried looking for papers about extending the PAM50, but could not find anything. I would greatly appreciate any advice.

# Answer
There is a CSV table in this paper with 33 sets of genes.
> 3 votes

---
Tags: gene, cancer, deep-learning
---


# Question
Title: Can ortholog genes have different names among different species?
To be specific, we have different lineages of E. coli. I wondered if the ortholog genes among them should all have the same name in all lineages or can they also have different names?

# Answer
Gene naming is pretty subjective and even a single gene can have many names. If you're lucky orthologs will be named the same in different lineages, but you can't rely on it. Better to use gene ids (e.g. from Ensembl) and an ortholog db to map the orthologs. However, if you're looking at different strains, this is more tricky since most db don't go below the species level
> 3 votes

---
Tags: phylogenetics, sequence-homology, orthologues
---


# Question
Title: How to make a video of a gromacs gmx dump standard output with VMD
I did `gmx trjconv` on GROMACS `trr` file to generate a `.xtc` file, and then used `gmx dump` and saved its standard output in a file. Then I tried to make a movie with VMD using the file. As a result, it does not work properly. The error in the terminal is as follows.

```
gromacsplugin) Cannot read header fromm '/Users/kondouakira/mount/Lab/1ns_sim/md_0_1_noPBC3.trj', no error

```

I use following command.

```
gmx trjconv -s md.tpr -f md.trr -o md_0.trr -pbc mol -ur compact -b 1 -force
gmx dump -f md_0.trr > md_0.trj

```

# Answer
> 2 votes
No need to use `gmx dump`, you can simply load a `xtc` or `trr` trajectory in VMD using:

```
vmd path/to/gro_file.gro other/path/to/trajectory_file.xtc

```

As always you can also use the GUI to first open your `gro` file, then right click on the molecule and select `Load data into molecule` before selecting your `xtc` or `trr` file.

After putting the representations as you see fit, if you want to export this as a movie (GIF, MP4, others...) you can go in Extensions -\> Visualization -\> Movie Maker. See Trajectories and Movie Making VMD official tutorial for more details.

# Answer
> 0 votes
I don't think that the output of `gmx dump` can be read by VMD: the goal of this command is to provide a human-readable output, not to convert one trajectory format into another.

`trj` files are **NOT** the output of `gmx dump`. `trj` files were used in "old" versions of GROMACS to store trajectories (see trj file format documentation of GROMACS 5), and this is what VMD is awaiting, not the output of `gmx dump`. This is also why it tells you the header is wrong.

In recent versions of GROMACS `trj` files are no longer supported, probably as they were not in a binary format, hence took too much disk space (see Trajectory files formats in GROMACS 2021 documentation). Converting a `.trr` in a `.trj` was done by `gmx trjconv`, but I don't think this is possible anymore (see `gmx trjconv` supported trajectory formats in GROMACS 5 vs GROMACS 2021).

Why would you absolutely need `gmx dump` output to be read by VMD? The only reason I can see is a handmade script doing modifications to your trajectory. In this case I would recommend to instead look into some molecular dynamics python libraries, like MDAnalysis, PyTraj or MDTraj, which will also allow you to make minute modifications to the trajectory file while still reading it in its original format. You can also look harder into GROMACS tools, as they are quite complete.

---
Tags: gromacs, molecular-dynamics
---


# Question
Title: getGeneLengthAndGCContent in EDASeq complains about "Ensemble gene ID"
I am trying to use the `getGeneLengthAndGCContent` function from the `EDASeq` library to retrieve gene lengths for c. elegans. I retrieved Ensemble gene ids from Biomart for c.elegans and I'm using those as input to getGeneLengthAndGCContent (I picked "ensemble\_gene\_id" from the filters of the worm annotation Mart). Here's an example of how I'm trying to use the function:

```
getGeneLengthAndGCContent("WBGene00001042", org, mode=c("biomart", "org.db"))

```

When I try to run this, I get

```
Error in getGeneLengthAndGCContent("WBGene00001042", org, mode = c("biomart",  :
  Only ENTREZ or ENSEMBL gene IDs are supported.

```

# Answer
> 0 votes
I found my answer here. It uses `getGeneLengthAndGCContent` but changes a small part of the function.

---
Tags: ensembl
---


title: Bioconductor Conda release delay

text:
# Question
Title: Bioconductor Conda release delay
How long does it take for the latest Bioconductor to end up on Conda? I think I saw a blog post about the challenges with this process, but I cannot find it. For example, Bioconductor 3.13 was released a week ago, but is still not available. For comparison, many CRAN packages tend to show up within a day.

Additionally, is there a way to keep track of that progress?

# Answer
> 4 votes
I started the roll-out on the 30th in the evening. There's not much of a delay for the October Bioconductor release, but for the summer release we first have to build the newest R version and then rebuild all conda CRAN packages. It's probably not widely known, but Bioconductor packages are only compatible with a single R version. For example, the ones being released now are not compatible with R 4.0 and will not be compatible with R 4.2, but will instead only ever work with R 4.1. This is regardless of what the metadata provided by Bioconductor package authors might say, the Bioconductor core team has requested this (there's actually not even a guarantee that packages built under R 4.1.0 will work with 4.1.1, but such issues are exceedingly rare). This is why there is a multi-week delay every summer.

On the plus side, at least we're not waiting for the first bug-fix release of R to start the migration anymore.

For Bioconda we've never gotten around to making a migration tracker. If you have a local clone of Bioconda, you can run `./scripts/bioconductor/rootNodes.py --sinceNDays 2 --rootNodes` to get a count of the packages completed in the past 2 days and a list of the root nodes on the dependency DAG that still need to be handled. That's how I estimate the fraction finished (currently about 600 packages).

For tracking the status on conda-forge you can see their status page: https://conda-forge.org/status/#r410

# Answer
> 3 votes
Devon will probably add a more precise answer, until then I just link his Tweet from last week:

https://twitter.com/dpryan79/status/1394714988720308226

Edit: 31.5.21

Rollout started yesterday, should take a week to complete:
https://twitter.com/dpryan79/status/1399273322601459712

---
Tags: bioconductor, conda, bioconda, anaconda
---


# Question
Title: rmsd between re-docked complex and co-crystallized complex
I would like to calculate the rmsd between re-docked complex and co-crystallized complex for docking validation. How do I calculate the rmsd using Pymol? I tried the following command

```
align re-docked_complex, co-crystallized_complex

```

But the above command gives the rmsd of two proteins. I would like to get the rmsd of re-docked and co-crystallized ligand. How can I achieve this using pymol?

# Answer
Welcome.

The Protein-Ligand-complex is the complete protein with Ligand!

If you only want to compare the two ligands, you can delete the rest of the proteins, or specify the input so only the ligand is choosen. (select Ligand and save selection. then use the function on that selections)

But be aware that for Docking comparison you probably don't want to align! but rather calculate the RMSD directly (because the absolute position and not the relative position is of interest?)

Maybe use: "rms\_cur" (https://pymolwiki.org/index.php/Rms\_cur)
> 1 votes

---
Tags: docking, pymol
---


title: Consensus sequence from vcf file

text:
# Question
Title: Consensus sequence from vcf file
I am generating a consensus sequence for SARS-Cov-2. I am visually inspecting every candidate SNP in a pileup.

Sometimes I find that there are two consistent variants for the same position, but I can only generate one consensus sequence.

Which variant should I target:

* The majority variant
* The minoritarian one if I find that this SNP is under-reported in other public samples?
* Are both required

What one is best practice?

# Answer
Making a consensus sequence tends to involve some arbitrary decisions. In this case, the common practices are:

1. Use the ancestral allele (this isn't always apparent)
2. Use the most common allele
3. Randomly choose

I think 1 and 2 are generally preferred, with 2 being easier to implement.
> 4 votes

---
Tags: vcf, variant-calling, covid-19, sars-cov-2
---


title: Which format is the most raw format

text:
# Question
Title: Which format is the most raw format
I wanted information regarding which format is considered the most raw from an illumina sequencer ?

1. fasta
2. fastq
3. bcl
4. bam
As per my research it should be **bcl** but I am not sure.

# Answer
BCL (binary base call) files are the most raw data that you will likely get. As a matter of fact, the images taken by the machine each cycle are even more raw but it makes no sense to interact with this data directly. The base calls are derived from counting the spots (clusters) on the images.

For any useful application you would convert them to fastq (e.g. using illumina `bcl2fastq`).

What is the goal of your question? Do you suspect some problems with your data?
Adding this information will make it easier to interpret what kind of answer you are looking for - and less likely that your question is downvoted.
> 1 votes

---
Tags: ngs, file-formats, illumina
---


# Question
Title: Batch effect correction using a subset of samples - using DE genes
I have some RNA-seq data with two very obvious batches as you can see in the PCA:

The samples of interest (A - H) are from tumor tissues. In addition, there are data for two cell lines I and J. The cell lines should be the same across the batches, while for the tumor samples there may be some additional biological difference.

My question is, if there is a way to use the cell line data to estimate the technical part of the batch information since there should very little biological difference here. Then I would like to use this to correct the technical aspect of the batch effect in the tumor samples.

When I run normal batch correction like `ComBat` or add the batch to the DE model the batch will be estimated over all samples including biological differences.

## Update

I figured I could determine which genes are differentially expressed between the batches in the cell lines only. As mentioned before, the cell lines are highly similar biologically. Thus these DE genes should represent the technical aspect of the difference.

The plot below shows the PCA when after removing genes DE in the cell lines. Obviously, the differences for the cell lines are greatly reduced. Also the tumor samples are now less clearly separated.

Most of the removed genes were from shorter RNA species and a look at sequencing QC later confirmed that in fact the fragment size in the second batch was lower. Therefore I'm positive that I could identify and remove the technical component of the batch effect.

## tl;dr

My data contains a strong difference between batches even for samples that are biologically the same (untreated cell lines). I looked at differentially expressed genes in the cell lines. This revealed to me that the second batch contained a lot more short RNAs due to lower fragment size of the library.

In the end I'm still wondering if there is an elegant solution for this requiring less footwork

# Answer
The batch effect is expected to be the same across all samples - therefore you want to estimate its effect across all of the samples. Just add batch as an additional parameter to the linear model (this is presuming samples sequenced in the two batches have been randomised).

You can show the effect of "batch effect" removal with a PCA plot by plotting the residual error, e.g. in R if using DESeq:

```
mypca <-  prcomp(t(assay(varianceStabilizingTransformation(dds)))
pc.resid <- resid(aov(mypca$x~batch.info))
ggplot(pc.resid,aes(x=PC1,y=PC2) + geom_point()

```

(I haven't checked the code, but it will be something like this).
> 1 votes

---
Tags: rna-seq, batch-effects
---


# Question
Title: Entrez (Biopython) esearch and efetch not returning sequence as expected
I'm trying to use Entrez (through Biopython) to download the sequence of a TMV replicase gene.

I have the following code:

```
# Lookup ID
search = Entrez.esearch(db='gene', term='Tobacco mosaic virus[Orgn] AND replicase')
read = Entrez.read(search)
idlist = read["IdList"]

# Fetch sequence
search = Entrez.efetch(db='gene', id=idlist[0], retmode='text', rettype='gb')
read = SeqIO.read(search, "genbank")
tmv_replicase = read.seq

```

But when I run it it throws the exception on the `SeqIO.read` line: `ValueError: No records found in handle`.

The IDs provided in idlist are: `['1494081', '1494082']`.

# Answer
Did you try to fetch from the database `nucleotide` instead of `gene` ?

```
# load modules
from Bio import SeqIO
from Bio import Entrez

# Lookup ID
search = Entrez.esearch(db='nucleotide', term='Tobacco mosaic virus[Orgn] AND replicase')
read = Entrez.read(search)
idlist = read["IdList"]

# Fetch sequence
search = Entrez.efetch(db='nucleotide', id=idlist[0], retmode='text', rettype='gb')
read = SeqIO.read(search, "genbank")
tmv_replicase = read.seq

```

It returns

```
>>> read
SeqRecord(seq=Seq('GTATTTTTACAACAATTACCAACAACAACAAACAACAAACAACATTACAATTAC...GGG'), id='MN912489.1', name='MN912489', description='Tobacco mosaic virus isolate TMV:BR:SC:02:01, complete genome', dbxrefs=[])
>>> tmv_replicase
Seq('GTATTTTTACAACAATTACCAACAACAACAAACAACAAACAACATTACAATTAC...GGG')

```
> 2 votes

---
Tags: python, covid-19, genbank, entrez, virus
---


# Question
Title: Which sequencing technologies are considered short read
I couldn't find any information on this online. Which of the following sequencing technologies are considered **short read** technologies? Illumina, Oxford Nanopore, Ion Torrent, Pacific Biosciences

# Answer
As mentioned by Greg, list of short-read sequencing platforms is listed on wikipedia.

Of those you mentioned Illumina and Ion Torrent are short-read sequencing platforms.

Usually what makes short reads short is a PCR step in the library prep or directly on the flowcell to amplify the signal, which is the limitation for how long stretches of DNA are possible to sequence. Long read sequencing technologies sequence directly isolated single molecules.
> 1 votes

---
Tags: ngs, homework
---