# Question
Title: What should GC coverage bias plot of exome data look like?
I want to see if there is any GC coverage bias in some paired-end Illumina human exome data. **Update** This is for a variant calling project. I used the `CollectMultipleMetrics` from GATK4 toolkit with `-PROGRAM` CollectGcBiasMetrics. An example of the command I used is:

```
gatk CollectMultipleMetrics -I input.bam -O prefix_dir_name -R input.ref_genome -PROGRAM CollectGcBiasMetrics -PROGRAM CollectAlignmentSummaryMetrics -PROGRAM CollectInsertSizeMetrics -PROGRAM QualityScoreDistribution -PROGRAM MeanQualityByCycle -PROGRAM CollectBaseDistributionByCycle

```

**Update** This is my result for the QC coverage bias with one sample highlighted:  There seems to be a large bias. I would expect the peak to be around 50% for human. And those weird peaks around 90% are also puzzling to me. On the other hand, I am working with exome data so it might be logical to have a bias towards the higher values of GC% because you have reads on GC-rich exons?

**Update** I have another plot with the bins produced by the above mentioned tool (.gc\_bias.pdf)

**Update** first part of useful file (gc\_bias.detail\_metrics) which shows the reads per bin, as suggested by @Maximilian Press. This I will plot.

```
## htsjdk.samtools.metrics.StringHeader
# CollectMultipleMetrics --INPUT WESP_run/results/aligned/RM8398_MF_S8_MD.sorted.bam --OUTPUT WESP_run/results/aligned/post_alignment_QC/multiple_metrics/RM8398_MF_S8 --PROGRAM CollectAlignmentSummaryMetrics --P
ROGRAM CollectBaseDistributionByCycle --PROGRAM CollectInsertSizeMetrics --PROGRAM MeanQualityByCycle --PROGRAM QualityScoreDistribution --PROGRAM CollectGcBiasMetrics --REFERENCE_SEQUENCE /mnt/scratch_dir/deute
koe/projects/HP/WESPipe/workflow/resources/reference_genomes/hg38/alt_plus/GCA_000001405.15_GRCh38_full_plus_hs38d1_analysis_set.fna --ASSUME_SORTED true --STOP_AFTER 0 --METRIC_ACCUMULATION_LEVEL ALL_READS --IN
CLUDE_UNPAIRED false --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_sec
rets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false
## htsjdk.samtools.metrics.StringHeader
# Started on: Sat Dec 03 14:25:15 CET 2022

## METRICS CLASS        picard.analysis.GcBiasDetailMetrics
ACCUMULATION_LEVEL      READS_USED      GC      WINDOWS READ_STARTS     MEAN_BASE_QUALITY       NORMALIZED_COVERAGE     ERROR_BAR_WIDTH SAMPLE  LIBRARY READ_GROUP
All Reads       ALL     0       139942  1326    19      0.413544        0.011357
All Reads       ALL     1       102412  61      17      0.025996        0.003328
All Reads       ALL     2       118957  74      16      0.02715 0.003156
All Reads       ALL     3       152181  929     20      0.266429        0.008741
All Reads       ALL     4       167732  248     16      0.06453 0.004098
All Reads       ALL     5       173319  232     16      0.058421        0.003836
All Reads       ALL     6       196365  304     15      0.067567        0.003875
All Reads       ALL     7       211217  343     17      0.070875        0.003827
All Reads       ALL     8       241125  441     18      0.079822        0.003801
All Reads       ALL     9       270272  1068    18      0.172463        0.005277

```

I have read a paper that uses the CollectGcBiasMetrics on exome data and they magically produce a plot without bias. I say magically because they don't explain how they use CollectGcBiasMetrics.

**Update** I added the target coverage plot.

**Update** Am I missing something? Should I, for instance, use my target regions in instead of the reference genome when using `CollectGcBiasMetrics`? I will try it, but I am not sure this is the correct thing to do... It is hard to find any workflow for exome data with gatk and particularly `CollectGcBiasMetrics`, which specifies to use the reference genome.

From the comment @Maximilian Press

> Do you have on-instrument diagnostics? Did you run FastQC on the raw reads? What do the other picard metrics look like? Did you manually inspect your alignments (strongly recommended)? Do you actually have coverage of the exome?

I will see if I can get some instrument diagnostics. I did look at the alignments, and as far as I can tell I don't see anything weird? Just that it is very deeply sequenced. I do see some intronic regions, with some reads, but only a few reads compared to the many on the targets, which I think is a consequence of the (too) deep sequencing. The fastqc of the raw reads didn't look troubling, and it did not give me any warnings. The picard metrics I thought did not look out of the ordinary. There is good target coverage. Yes, it is sequenced quite deep. Also, the lab tech did tell me they over-cycled because they wanted to be sure to have all the regions/targets.

**Update** So this bias might be due to too deep sequencing? And another important question, can I still reliably use this data? If so, I guess I can perform some corrections like suggested here

Any help/suggestions would be very welcome. Thank you in advance.

---

# Answer
> 2 votes
> Am I missing something? Should I, for instance, use my target regions instead of the reference genome?

Yes you should do this. Why do you think it's wrong doing that? It should be the first thing to check. And that GC bias might be effected by multiple factors. Most prominent one would be PCR over amplification. Most PCR polymerases give horrible GC bias when library is over cycled. There are many articles published on this.

I also suggest that you look at the number of reads in each bin. Peaks around 90% GC should not puzzle you so much. There is chance that you have only few reads in these bins out of millions of reads.

# Answer
> 1 votes
Just to formalise my comments.

@Supertech's is the right answer.

**Background** The one thing that is unusual is that the bias peaks at around 75%. The only time I've seen strong GC% bias is for particular groups of bacteria. 75% is about the maximum bias that can be achieved in a *protein coding gene* and these "extreme bacteria" hit that limit.

I better try and explain that a bit better. The primary selection pressure on a gene is the amino acid composition and the GC% is a second issue, i.e. selection pressure, for some organisms (e.g. certain bacteria, insects ..). Thus amino acids are not mutating, but the triplet codon degeneracy permits a bias in the GC% without alternating a given amino acids. In other words this GC% bias is occurring only for "neutral" mutations (clearly there must be a selection pressure, so they are not really neutral).

**Summary and personal thoughts** In summary the GC% limit is around 75% for a protein coding gene and beyond that it requires the amino acids to be changed. The peak in @Dandelion's data on that limit is weird because GC% artefacts from *Taq* errors are not subject to selection pressure: there should not be a peak.

"Contamination" I would avoid as an explanation because in my experience these groups of bacteria are unlikely to invade a lab environment. Sample mix up is a possible explanation. However, I would be cautious about assigning blame because in the chain of command within the wet-lab it could have occurred anywhere.

**Official response** To reiterate the formal and official answer would be @supertech's response and the 75% GC peak remains a mystery.

---

**Please note** The second plot has a distribution around GC content around 40% ... There appears a discrepancy between plots. Neither output appears compatible with the GC% of humans, but 40% GC with a reasonable proportion of that data on the 50% mark is better the 75% peak plot. Personally I would investigate the empirical basis for the differences.

---
Tags: gatk, illumina, human-genome, coverage, exome
---


# Question
Title: Unable to open .bam file in C++ using SeqAn due to 'seqan::UnknownExtensionError'
I am trying to open .bam files in C++ to extract reads occurring at specific scaffolds and loci. I essentially want to call "samtools view sample.bam -o sample.sam scaffold:pos-pos" from C++. I have tried system("samtools view sample.bam -o sample.sam scaffold:pos-pos") and can explain the issues if you want to know, but for the sake of conciseness I will leave the details out of this post. So, I am using SeqAn. When I run my script, I get the following error:

```
terminate called after throwing an instance of 'seqan::UnknownExtensionError'
  what():  Unknown file extension of copyPolly_SRR6511930.bam: iostream error

```

I am running my program on my institution's supercomputer clusters (Texas A&M Terra cluster), to which I do **not** have sudo access. I have loaded the following modules:

```
GCC/9.3.0
SeqAn/2.4.0-GCCcore-9.3.0
zlib/1.2.11-GCCcore-9.3.0

```

Here is my minimal script:

```
#include <seqan/bam_io.h>
#include <zlib.h>

using namespace seqan;

int main(){
    CharString bamFileName = "copyPolly_SRR6511930.bam";
// Open input file.
    BamFileIn bamFileIn(toCString(bamFileName));
    if (!open(bamFileIn, toCString(bamFileName))){
        std::cerr << "ERROR: Could not open " << bamFileName << std::endl;
        return 1;
    }

    unsigned numUnmappedReads = 0;

    try{
        // Read header.
        BamHeader header;
        readHeader(header, bamFileIn);

        // Read records.
        BamAlignmentRecord record;
        while (!atEnd(bamFileIn)){
            readRecord(record, bamFileIn);
            if (hasFlagUnmapped(record))
                numUnmappedReads += 1;
        }
    }
    catch (Exception const & e)
    {
        std::cout << "ERROR: " << e.what() << std::endl;
        return 1;
    }
    return 0;
}

```

I have tried the following solutions:

-Adding `#define SEQAN_HAS_ZLIB (1)` to my script's header

-Adding `-D SEQAN_HAS_ZLIB` to my compilation such that it reads `g++ script.cpp -D SEQAN_HAS_ZLIB`

-Adding `-lz` to my compilation such that it reads `g++ script.cpp -lz`

I am not dead-set on using SeqAn for this project, so if anyone knows of another C++ library which extracts specific reads from .bam files, please share.

Update: I've never used HTSlib and it seems to have a fairly steep learning curve, so if anyone has an answer for the SeqAn error, that would still be very much appreciated.

# Answer
> 3 votes
The helpdesk for my institution's supercomputers showed me the following flags:

```
g++ seqan_test.cpp -DSEQAN_HAS_ZLIB -lz -lpthread

```

This fixed the error. I'm now running into different errors and am suspicious that I am getting them because I am using SeqAn2 rather than SeqAn3. I'm going to try HTSlib.

# Answer
> 3 votes
Htslib is the official library to access alignment files. It supports more formats. The typical way to read through a SAM/BAM/CRAM:

```
#include "htslib/sam.h"
int main(int argc, char *argv[])
{
  samFile *fp = sam_open(argv[1], "r");
  bam_hdr_t *h = sam_hdr_read(fp);
  bam1_t *b = bam_init1();
  while (sam_read1(fp, h, b) >= 0) {
    if (b->core.tid < 0) continue;
    int end_pos = bam_endpos(b);
    printf("%s\t%d\t%d\n", h->target_name[b->core.tid], b->core.pos, end_pos);
  }
  bam_destroy1(b);
  bam_hdr_destroy(h);
  sam_close(fp);
}

```

PS: the above started as a comment for "relevant but minor or transient information to" the original question. I didn't consider htslib a direct answer to the question, but because this had been changed to an answer, I added some details anyway. By the way, seqan2 has largely been discontinued. Use seqan3.

# Answer
> 1 votes
I'm not quite sure why assigning a string to a variable would lead to a complaint about file extensions. From the SeqAn BAM tutorial, it looks like you need to use the *BamFileIn* object for reading BAM/SAM files:

```
#include <seqan/bam_io.h>
using namespace seqan;

int main(){
    CharString bamFileName = "copyPolly_SRR6511930.bam";
    // Open input file, BamFileIn can read SAM and BAM files.
    BamFileIn bamFileIn;
    if (!open(bamFileIn, toCString(bamFileName))) {
        std::cerr << "ERROR: Could not open " << bamFileName << std::endl;
        return 1;
    }
    return 0;
}

```

Edit: after seeing the updated code... I've got nothing, sorry. It all seems like reasonable code that shouldn't be generating errors.

# Answer
> 1 votes
I think `SEQAN_HAS_ZLIB` flag not works correctly as it supposed to be. So you can check the seqan library to see if `SEQAN_HAS_ZLIB` works correctly, and seek the reason of extension error.

Seqan library is a mature library which using template to define general-purpose classes and structs.

1.The reason for extension error.

Step into the seqan library:

When you call BamfileIn: `BamFileIn bamFileIn(toCString(bamFileName));`

An object of BamFileIn struct is created, and `bamFileName` file is opened by the struct construction function.

And when this object is created, some extra information is wrapped with the object to check the type of file.

This information include suffix categories. In bam\_file.h header file.

```
template <typename TDirection, typename TSpec>
struct FileFormat<FormattedFile<Bam, TDirection, TSpec> >
{
#if SEQAN_HAS_ZLIB
    typedef TagSelector<
                TagList<Bam,
                TagList<Sam
                > >
            > Type;
#else
    typedef Sam Type;
#endif
};

```

As you can see, if `SEQAN_HAS_ZLIB` is setted, "Bam" file will be added to TagSelector which will be wrapped into the BamFileIn struct and redefined by `typedef` as `format`.

Then \_open() function is called to open Bam file, and extension is checked in guess\_format.h

```
template <typename TFilename, typename TTagList>
inline bool guessFormatFromFilename(TFilename const &fname, TagSelector<TTagList> &format)
{
    typedef typename TTagList::Type TFormat;

    if (value(format) == -1 || isEqual(format, TFormat()))
    {
        // if tagId is set to -1 (auto-detect) or the current format (TFormatTag) then test for TFormatTag format
        if (guessFormatFromFilename(fname, TFormat()))
        {
            assign(format, TFormat());
            return true;
        }
    }
    return guessFormatFromFilename(fname, static_cast<typename TagSelector<TTagList>::Base &>(format));
}

template <typename TFilename, typename TFormat_>
inline bool guessFormatFromFilename(TFilename const & fileName, Tag<TFormat_> /*formatTag*/)
{
    typedef typename Value<TFilename>::Type                                     TValue;
    typedef ModifiedString<TFilename const, ModView<FunctorLowcase<TValue> > >    TLowcase;
    typedef Tag<TFormat_>                                                       TFormat;

    TLowcase lowcaseFileName(fileName);
    for (unsigned i = 0; i < sizeof(FileExtensions<TFormat>::VALUE) / sizeof(char*); ++i)
        if (endsWith(lowcaseFileName, lowerString(FileExtensions<TFormat>::VALUE[i])))
            return true;

    return false;
}

```

The above code is two instance of `guessFormatFromFilename` function. It's pretty complex here, but if `SEQAN_HAS_ZLIB` is set correctly, `TagSelector` will include bam and sam, and the first `guessFormatFromFilename` above will return true in if condition, and the error will not happend. But if `SEQAN_HAS_ZLIB` is not working, second formt `guessFormatFromFilename` function will be called and only "sam" suffix is exist in `TFormat`, and it will be exhibited by `FileExtensions<TFormat>::VALUE[i]`. You can using std::cerr to print the suffix on your screen to see if there is only "sam" suffix is exist in `TFormat`. Althrough you add -DSEAN\_HAS\_ZLIB, I think you could check the code in library to see if it really works fine.

If `SEQAN_HAS_ZLIB` is working, you could check other variables or functiones related with open bam process.These function are in `bam_file.h`, `formatted_file.h`, `fundamental_tag.h` and so on.

---
Tags: sam, samtools, c++
---


# Question
Title: ERR: error while loading shared libraries: libjulia.so.1: cannot open shared object file [Julia] - Ubuntu 22.04
I'm trying to build a Dockerfile. It has a tool called Atria that utilizes Julia. Since Ubuntu 22.04 does not have Julia package, I had to resort to installing it. The following is the Dockerfile:

```
FROM ubuntu:22.04

# Install tzdata
RUN apt-get update &&\
    DEBIAN_FRONTEND=noninteractive TZ=Etc/UTC apt-get -y install tzdata

# Install Python, Bowtie2 and Java
RUN apt-get install -y python3.10 python3-pip \
    openjdk-8-jdk  \
    bowtie2 \
    wget

RUN  apt-get install --yes --no-install-recommends \
     zlib1g-dev \
     libbz2-dev \
     liblzma-dev

# Install RSeQC
RUN pip3 install RSeQC

# Install biopython=1.80
RUN pip3 install biopython

#Install Julia
RUN wget https://julialang-s3.julialang.org/bin/linux/x64/1.8/julia-1.8.1-linux-x86_64.tar.gz && \
    tar zxvf julia-1.8.1-linux-x86_64.tar.gz && \
    mv julia-1.8.1/bin/julia /usr/local/bin/julia && \
    chmod +x /usr/local/bin/julia

RUN /sbin/ldconfig -v

# Install Atria
RUN wget https://github.com/cihga39871/Atria/releases/download/v3.1.2/atria-3.1.2-linux.tar.gz && \
    tar -zxf atria-3.1.2-linux.tar.gz && \
    mv atria-3.1.2/bin/atria /usr/local/bin/atria && \
    chmod +x /usr/local/bin/atria

#Atria dependencies
RUN  apt-get install pigz pbzip2

# Install findtail
 RUN wget https://storage.googleapis.com/google-code-archive-downloads/v2/code.google.com/findtail/findtail_v1.01 && \
     mv findtail_v1.01 /usr/local/bin/findtail_v1.01 && \
     chmod +x /usr/local/bin/findtail_v1.01
# Cleanup
RUN apt clean

```

What I have tried so far is :

1. Julia v1.6.7 and v1.8.4
2. `RUN /sbin/ldconfig -v` as suggested here

How to fix it? Should I try to install Julia using PPA?

edit- The error:

```
root@b0f4b17a7c18:/# atria
/usr/local/bin/julia: error while loading shared libraries: libjulia.so.1: cannot open shared object file: No such file or directory

```

It is recognizing the atria command tho.

# Answer
> 2 votes
You just need to make sure you include the required library files when you install Atria. Since you're not building any software (i.e. using `make`), the usual way is to simply move the directory under `/opt` and add the executable(s) to your `PATH` using the `ENV` keyword. Note that Julia comes bundled with Atria, so all you need is:

```
FROM ubuntu:22.04

RUN apt-get update && apt-get install -y \
    bowtie2 \
    openjdk-8-jdk \
    pbzip2 \
    pigz \
    python3.10 \
    python3-pip \
    python3-biopython \
    wget \
    && apt-get clean \
    && rm -rf /var/lib/apt/lists/*

# Install RSeQC
RUN pip3 install --no-cache-dir RSeQC

# Install Atria
RUN wget -q https://github.com/cihga39871/Atria/releases/download/v3.1.2/atria-3.1.2-linux.tar.gz \
    && tar xf atria-3.1.2-linux.tar.gz \
    && mv atria-3.1.2 /opt/atria

# Install findtail
RUN wget -q https://storage.googleapis.com/google-code-archive-downloads/v2/code.google.com/findtail/findtail_v1.01 \
    && chmod +x findtail_v1.01 \
    && mv findtail_v1.01 /usr/local/bin/findtail

# Add the Atria executable to the PATH
ENV PATH="/opt/atria/bin:${PATH}"

```

# Answer
> 1 votes
Just to note `RUN apt clean` ... is much better performed using `&& apt clean`, like in @Steve's approach, recommended for complicated reasons.

**Suggestion**: what about using the official Julia Docker as the base rather than create it from scratch? There isn't a Ubuntu Julia but there is a Debian base.

```
docker run -it --rm Julia

```

or

```
docker pull julia

```

or via a Dockerfile,

```
FROM Julia:buster

```

or

```
FROM Julia:latest

```

---
Tags: docker
---


# Question
Title: What is the correct method of identifying the interacting residues for specific molecular docking?
I'm trying to do a specific molecular docking using Autodock Vina in PyRx. So, one way to do that is to mark all of the interacting residues and make sure the grid box encompasses all of that interacting residues.

So, my protein already has a co-crystalized ligand attached to it. However, I'm not sure what method to use.

After searching online, I found out that I should identify the interacting residues of that ligand. Then, I have opened this protein-ligand complex in 2 applications and 1 web server.

Here are the information of the protein-ligand complex:

* The protein PDB ID is: 5TL8 (Naegleria fowleri CYP51-posaconazole complex)
* The ligand is Posaconazole
* The cofactor and other small molecules: HEME, 1,2-ETHANEDIOL, calcium ion

I tried opening the complex in Discovery Studio Ligplus

and Protein Plus webserver (this server has more interacting residues). All of them have different residues.

Which one should I follow? Should we included the HEME or do we just look at the Posaconazole itself?

# Answer
> 3 votes
I think I've found the answer to my own question.

So, I tried using this PRANK web server for predicting the binding site of my protein. The results is suprisingly very similar to the one from the Protein Plus web server. So, I think it's better to use the one from Protein Plus.

You can see the results on the link below.

References:

1. https://prankweb.cz/viewer?id=5tl8&database=v3-conservation-hmm&created=2022-01-27T07:50:59&lastChange=2022-01-27T07:50:59%20%20&structureName=structure.pdb&predictionName=5TL8&predictedStructure=undefined
2. https://proteins.plus/5tl8

# Answer
> 4 votes
Not an answer either. But a series of caveats.

First, `HEM` parameterisation. For molecular thermodynamics tools you will need to generate a topology file for this bad-boi. The iron has 2 donor bonds and 4 single bonds, one of which with the coordinating cysteine and the other with an oxygen, which will attach the substrate, in this case a steroid —let's ignore the radical state of the activated bound oxygen, but do keep in mind that the partial charges and atom type auto-assigned will be wrong. A cheat would be constraining said molecular oxygen and using a porphyrin ring with no iron: not reality, but would make the experiment work, is as physically questionable as any iron-bound topology you'll get and the oxygen is the key player here.

Second, cytochromes act on greasy substrates, there are few if any polar interactions generally, so constraining to polymer residues may not be that useful. Specificity is controlled by the gating and entranceway (think of trucks reversing into a loading bay), which is annoying the gating loops are floppy loops that are generally not crystallised, so tethered MD from a bound to unbound state (undocking) is sometimes done. There's a lot of literature on P450 and P411 cytochromes, but often it's not as clear cut and simple as say a hydrolase.

However, if your substrate you want placed is similar to one present you may be interested in a tool of mine Fragmenstein, which can be used to convert a compound to a similar one (stripping the haem and adding it later: I had a Rosetta params file for `HEM` somewhere but not sure where it is nowadays)

# Answer
> 3 votes
not an answer, just more inputs, but https://proteins.plus/5tl8#poseview

gives me this:

references:

ProteinsPlus: interactive analysis of protein–ligand binding interfaces

Molecular complexes at a glance: automated generation of two-dimensional complex diagrams

Drawing the PDB: Protein−Ligand Complexes in Two Dimensions

from this last paper :

> Hydrogen bonds are implemented mainly following the measures published by Desiraju and Steiner in 2001.14 The optimal distance between two atoms connected by a hydrogen bond is set to 1.9 Å with a tolerance of 0.5 Å. Additional to this measure, the acceptor−hydrogen−donor angle must not fall below 120°. Hydrogen atoms that are bound to a noncarbon atom are treated as hydrogen donor candidates. Potential acceptor atoms are either nitrogen, oxygen, or sulfur atoms provided that they are uncharged or negatively charged and that their surface is accessible.

> Metal interactions are calculated between metal atoms embedded in the receptor and metal acceptor atoms, which are identical to the hydrogen bond acceptor atoms. Their geometry is based on the calculated coordination geometry of the metal.15 Each coordination point that is not occupied by a receptor atom is checked for close ligand atoms, and the maximal distance deviation is set to 0.8 Å. If no geometry can be calculated, a sphere with a radius of 2 Å is placed around the metal. In this case, all atoms lying on the sphere, again with a tolerance of 0.8 Å, are regarded as potential interaction partners.

> In contrast to the formerly mentioned interactions, hydrophobic contacts are estimated based on the distance between two hydrophobic atoms only. They are visualized not by a dashed interaction line but by drawing the label of the contacting residue and a spline segment denoting the hydrophobic part of the ligand. Because many atoms typically form a hydrophobic subpocket, this representation reflects the interaction geometry better. A prerequisite for a hydrophobic contact is that at least three hydrophobic ligand atoms lie in the range of the currently examined receptor residue. Hydrophobic atoms are in this context carbon atoms with accessible surface and halogens. The maximum distance is set to the sum of the van der Waals radii of atoms in question and a tolerance of 0.8 Å.

**EDITED :**

Adding more on a side note I came accross this paper, I am completely naive about the topic but the title :

> “Observing Noncovalent Interactions in Experimental Electron Density for Macromolecular Systems: A Novel Perspective for Protein-Ligand Interaction Research.” Journal of chemical information and modeling vol. 62,7 (2022): 1734-1743. doi:10.1021/acs.jcim.1c01406

is kind of intriguing, I am not an expert and cannot vauch for its strictness but to me sounds like a good why to validate an interaction like seeing it even if in abstract terms:

> ... Non-covalent interactions (NCI) can be observed by pinpointing the ED saddle point, i.e. (3,−1) critical points and further quantified by measuring ED deviation from a homogeneous electron distribution using density and its first derivative (s = \[1/(2(3π2)1/3)\]|∇ρ|/ρ4/3)....

Please feel free to add to this, its not my field.

https://www.biorxiv.org/content/10.1101/2022.01.24.468575v2.full

---
Tags: protein-structure, pdb, docking, structural-biology
---


# Question
Title: CONECT doesnt include all the connections the 1vfl macro molecule have in a pdb file
I have a 1vfl macro molecule with 3000 atoms in it. When i open it in DC accelrys visualizer. It has much bonds. but when i open the pdb file it just have 8 connections. You can download the pdb here(I'm linking because the pdb is much) : https://www.rcsb.org/structure/1VFL
I'm using this pdb for my program so i need to know all the connections between atoms.

Any help would be good; I'm a complete beginner at using PDB.

# Answer
> 5 votes
The `CONECT` entries specify the bonding of the ligands (`HETATM`), not your polymeric residues (`ATOM`). By residue I mean a collection of atoms that forms a unit, not only your standard amino acid definition.

You can export in PyMOL all connections if you wanted, but the connectivity of the polymeric residues is standardised and does not need to be specified as they ought to obey the definitions (the 'chemical component' dataset), for example https://www.rcsb.org/ligand/ILE. If I were to have a ligand, `LIG`, it would be `HETATM` entered and with `CONECT` entries as most likely it is not the RCSB's https://www.rcsb.org/ligand/LIG.

Crosslinks between certain residues ought to be entered via the `LINK` record, not `CONECT` so they will not feature.

Some programs connect atoms by proximity regardless of stated connectivity —spiderwebs in PyMOL.

An odd corner case are non-standard amino acids, such as phosphorylations (e.g. https://www.rcsb.org/ligand/SEP) which are `ATOM` not `HETATM`, which are not meant to have `CONECT` entries, but some were deposited with them. Water (`HOH`) (from neutron diffraction structures) does not need connections.

---
Tags: pdb
---


# Question
Title: Best method to compare differentially expressed genes between 2 single cell clusters for GSEA
I want to get a ranked gene list for GSEA analysis from single cell data, which needs a group differential expression statistic for ranking.

Which is the best method to compare the differently expressed genes between 2 clusters `scanpy.tl.rank_genes_groups` (https://scanpy.readthedocs.io/en/stable/generated/scanpy.tl.rank\_genes\_groups.html) or `MAST.zlm` (https://rdrr.io/bioc/MAST/man/zlm.html) or `FindMarkers` function of `Seurat` (https://satijalab.org/seurat/reference/findmarkers)?

Thanks

# Answer
I wish to focus my contribution singly on ScanPy vs. Seurat FindMarkers.

ScanPy's claim is it is essentially a speeded up version of Seurat FindMarkers with better performance (discussed below) written in Python.

What they are doing are essentially datamining the expression signal using multivariate statistics (PCA) focused through tSNE. There is no controversy or dispute with this approach at all (please see **Caveat** I shouldn't use the term "datamining" in context - bad habit).

What ScanPy is doing is using Graph theory `Louvian groups` to extract the differential signal, presumably from the tSNE focused PCA. Cool.

If you've already got to grips with ScanPy then leveraging it as a data mining approach - to me - looks sensible. The statistics will not be hugely powerful, that doesn't mean they are bad, it approximately means they are solid.

**What is ScanPy doing?** The speeding up is via parallelisation, which is very easy in Python. The implication is Seurat is not parallelised, which would amaze me. The authors talk about some form of 'chunking' (they say 'partially loading the object') via HDF5. This is officially what they say

> ANNDATA is similar to R’s EXPRESSIONSET \[26\], but supports sparse data and allows HDF5-based backing of ANNDATA objects on disk, a format independent of platform, framework, and language. This allows operating on an ANNDATA object without fully loading it into memory

I am confused by this, because `HDF5` is a rapid I/O format - its not a datachunking format which his what the authors imply. Yes it facilitates datachunking but simply via rapid I/O. `HDF5` is of course independent Python, but why discuss that at length (we all know that) when the critical information is how the partially loaded object is achieved. There's a lot of strategies to achieve the last bit - none of which are discussed.

The overall levels of high performance Python (outside ScanPy) - because the reality is Python's core code is slow - are three steps:

1. Parallelisation and code optimisation.
2. Data chunking - using algorithm strategies to partial load data e.g. Hadoop strategies.
3. Compiling via `JIT`, `numba` and static typing via `Cython`.

The authors have done point 1, however again `HDF5` doesn't facilitate point 2 above, all it does is enable faster read/write. `HDF5` is not actually the fastest I/O either both `Pacquet` (sp?) and `feather` are faster but normally for `pandas` dataframes. `HDF5` is actually just as fast as `pickle` under many tests (!).

So the paper is a bit weird in that it is sparsely described - but this does often happen admittedly. The partial loading stuff via `HDF5` is weird, because they are not describing how that is done.

The other bit that is weird is the ML and ANNDATA class. They state,

> ANNDATA is similar to R’s EXPRESSIONSET \[26\]

I rarely use `R` so to me that isn't clear. The other weird thing is the authors talk about ML, they use an ML library, but give no further details nor why it is used. What I would assume is either: that they use it to assist with `pseudotemporal ordering` analysis **OR** they use it for unsupervised learning (PCA, tSNE). The latter is more likely. `Louvian groups` are from the `igraph` package they use

So in my personal view, despite the paper being a bit (to very) obscure I see no reason why ScanPy cannot replace Seurat purely in terms of the datamining being performed. Python is just as good at PCA and tSNE etc ... The use of Louvain groups is cool to trap that information. The key benefit of ScanPy is:

* your already familiar with it.
* faster - it definitely could be accelerated further.
* good for Python

The drawbacks I can see are:

* RNA analysis is dominated by `R` and you'll likely pick up a `R` focused reviewer (which could be fun).
* As *per* your previous post there appear shortcomings with ScanPy - but like the development team can't think of everything first time round.
* ScanPy2 is not released.
* its not generating a critical (test) probability to distinguish between expression levels just like Seurat - there are real advantages to this however.

**Caveat** I keep saying "data-mining" and I shouldn't. If the reviewers are perfectly happy with the above analysis alone, that's cool. This is because what I repeatedly call "data-mining" are other investigators final results, particularly in RNA expression. If your aim is to publish a big RNA-seq data set and identify the core features of e.g. a cellular immune response - good approach. I don't generate data, mine is purely analytics and algorithms, so my perspective (and reviewing standards) are different.

---

Just to state my interest in DEGs (or DGE). I use unsupervised learning, PCA, tSNE to generate a hypothesis and get a feel for the data. Then however, I do customised full-scale ML (supervised learning) to interrogate the data. Finally, I use DEG to assess which way around the ML result is pointing (under- or over-represented). I will not let the DEG override the ML (supervised learning) analysis - which I see as much stronger analysis. My context is different however and I'm consolidating the result is via extensive custom databases, i.e. very different.
> 1 votes

---
Tags: single-cell
---


# Question
Title: How to preserve colour assignments in matplotlib when a category is empty?
Some of my samples have empty clusters. I want to keep the same color for each cluster in different plots.

How do I get `sc.pl.umap` to keep the same colors for each cluster, even when a cluster is empty?

```
from matplotlib.pyplot import rc_context
with rc_context({'figure.figsize': (3, 3)}):
    fig, (ax1, ax2) = plt.subplots(2, 2, figsize=(6,6))
    ax1[0] = sc.pl.umap(adata[adata.obs['sample'] == adata.obs['sample'].cat.categories[0]],
                 color = 'leiden_r1', size=20,
                 legend_loc = 'on data', ax=ax1[0],
                 palette = mycolormap_26,
                 frameon=True, legend_fontsize = 8,title = adata.obs['sample'].cat.categories[0], show=False)
    ax1[0].set_xlabel('')
    ax1[0].set(xlim=xlim, ylim=ylim)
    ax1[1] = sc.pl.umap(adata[adata.obs['sample'] == adata.obs['sample'].cat.categories[1]],
                 color = 'leiden_r1', size=20,
                 legend_loc = 'on data', ax=ax1[1],
                 palette = mycolormap_26,
                 frameon=True, legend_fontsize = 8,title = adata.obs['sample'].cat.categories[1], show=False)
    ax1[1].set_xlabel('')
    ax1[1].set_ylabel('')
    ax1[1].set(xlim=xlim, ylim=ylim)
    ax2[0] = sc.pl.umap(adata[adata.obs['sample'] == adata.obs['sample'].cat.categories[2]],
                 color = 'leiden_r1', size=20,
                 legend_loc = 'on data', ax=ax2[0],
                 palette = mycolormap_26,
                 frameon=True, legend_fontsize = 8,title = adata.obs['sample'].cat.categories[2], show=False)
    ax2[0].set(xlim=xlim, ylim=ylim)
    #ax2[0].set_xlabel('')
    ax2[1] = sc.pl.umap(adata[adata.obs['sample'] == adata.obs['sample'].cat.categories[3]],
                 color = 'leiden_r1', size=20,
                 legend_loc = 'on data', ax=ax2[1],
                 palette = mycolormap_26,
                 frameon=True, legend_fontsize = 8,title = adata.obs['sample'].cat.categories[3], show=False)
    ax2[1].set(xlim=xlim, ylim=ylim)
    #ax2[1].set_xlabel('')
    ax2[1].set_ylabel('')

```

Thanks

# Answer
My recommendation is to report the issue directly on the `ScanPy` Github page here: https://github.com/scverse/scanpy/issues

This is a shortcoming with the testing of the algorithm, without any question. Replicates with experimental modification are standard procedure in this line of analysis. For example, one *in vitro* cell culture does not receive the drug therapy, i.e. a control, and the rest do, which results in one group with a complete absence of a transcript profile. The authors simply have not described this possibility in any of their documentation and therefore may not have originally considered it.

What I suspect has happened is the authorship only tested extensive and complete data - like really big data sets. "BigData" is what `ScanPy` is targeting because they look to have adapted to large data in Wolf *et al* Genome Biology 2018 paper, notably by specifying how quick their algorithm performs with increasing data set size.

What I suspect has happened is that smaller data sets, with experimental replication and modification, were not tested so these issues simply get missed. It's easy to overlook, the authors can't think of everything.

The authors of `ScanPy` are responding to issues arising unevenly. So it will be hit and miss whether they respond. Similar issues should have been reported to them before because its a common wet-lab experimental design. However, I've not checked their full intray.

They should provide the bit of code that will prevent the tedium described above, but then again they may not. If they do, please do write back. The whole `ScanPy` project is interesting.
> 1 votes

---
Tags: single-cell
---


title: Problems with BED format for ClinCNV

text:
# Question
Title: Problems with BED format for ClinCNV
I have some problems with ClinCNV(https://github.com/imgag/ClinCNV):

I made a BED file from a BAM, then i made a GC count with ngs-bits (link: https://github.com/imgag/ngs-bits ) The main idea of redacting a BED file is making GC content, but i think i need to remove all other columns. Here an example:

My BED:

```
chr1    10041   10086   NB501410:34:H23TGAFXY:1:11102:23474:15441/1     0       +       0.4889

```

In description says:

```
chr1    12171   12245   0.4595

```

I tried to run script with my BED file and there is an issue:

```
 Rscript ClinCNV/clinCNV.R --normal myfolder/myfile.cov --out myfolder/ --bed myfolder/myfileGCannotated.bed --folderWithScript ClinCNV/

```

and there is an output:

```
"We run script located in folder /root/myfolder . Please, specify ABSOLUTE paths, relative paths do not work for every machine. If everything crashes, please, check the correctness of this path first."

[1] "START cluster allocation."

[1] "Cluster allocated."

[1] "END cluster allocation."

[1] "We are started with reading the coverage files and bed files 2023-01-02 13:49:43"

Error in fread(fileName, header = F, stringsAsFactors = F, sep = "\t") :
  File 'myfolder/myfileGCannotated.bed' does not exist or is non-readable. getwd()=='/root/ClinCNV'
Calls: ReadFileFast -> as.data.frame -> fread
Execution halted

```

Now I started the script with absolute paths and it works but I have a new problem:

```
Error in setnames(localDf <- as.data.frame(fread(fileName, skip = 1, header = F,  :
  Can't assign 4 names to a 7 column data.table
Calls: ReadFileFast -> setnames
Execution halted

```

# Answer
The error message is telling you:

> "We run script located in folder /root/myfolder . Please, specify ABSOLUTE paths, relative paths do not work for every machine. If everything crashes, please, check the correctness of this path first."

So it sounds like you simply need to use absolute instead of relative paths. Something like this:

```
Rscript /home/yourUser/ClinCNV/clinCNV.R \
    --normal /home/yourUser/myfolder/myfile.cov \
    --out /home/yourUser/myfolder/ \
    --bed /home/yourUser/myfolder/myfileGCannotated.bed \
    --folderWithScript /home/yourUser/ClinCNV/

```

Also, don't run things as root unless absolutely necessary. if anything goes wrong or if any of your scripts or tools has a bug, you could end up breaking your machine if you run as root.
> 1 votes

---
Tags: ngs, bash, bed, coverage
---


# Question
Title: Should be RT-qPCR values standardized before PCA analysis?
I come from this question of gene exression analysis: Should PCA be standardized for gene expression?

My experiment is based in 151 x 51 (individuals/samples x genes), in which, patients are subjected to 3 possible groups (~ balanced groups), so we have 2 measurements of gene expression: one in the baseline, and one after the intervention. The gene expression is actually summarised in an scalar (number) called fold-change. Fold-change is calculated through the ***Ct*** (*cycle threshold*), which is the number of the cycle in which you are able to detect a certain amount of the target gen. This value or cycle is normalised using Ct values of 3 reference genes which are supposed to remain stable, and each patient is measured using oneself baseline value, in the same run, to remove run-to-run variation. As you can see I am working with RT-qPCR and relative quantification values.

As you can see, "normalization" or "standardization" is performed, not in the strict meaning of mathematical foundation, but to normalize quantitation of targets in the experiments. I was wondering if, taking into account that normalization is done, when I try to perform PCA to see how these genes behave, if they form cluster, if the intervention received has influence of the expression, patterns...

Reading about, I found quite divergent opinions about scaling or standardizing the results. The conflict, as far as I can see in my own experiment, lies in all my variables are genes, measured with the same high-throughput analyzer, different batches, but normalized with reference genes, and their own baseline values.

On one hand I found quite homogeneous working with scaled values, but on the other hand, I am working with values that seem to be working on the same scale, and have been through several filters to normalize results. Not sure if executing more mathematical operations could modify in the wrong way. **What is your opinion about it?** (all is about using the covariance matrix or correlation matrix based PCA, as far as I understand)

Furthermore, the dataset contain missing values, so I am working with the missMDA package to impute them. **In case I scaled the data, what should I do first: scale or impute the missing values?**

The steps to calculate fold-change:

**Update 1:** About the technology used, is a OpenArray plate (analyzed with QuantStudio™ 12K Flex Real-Tim) which contains subarrays with customised target genes.It is a 3,072-reaction high-throughput real-time PCR. The reaction are individually performed (as far as I know, getting *Ct* for each sample and assay). About the 2nd point, it is 2x time series study (baseline and post-intervention)- But if you see the equations to calculate the fold-change, I am going to obtain a single value of gene expression combining the 2 *Ct* or measurements for each patient (as I told: baseline and post-intervention).

About the "genes aim" you can see that my interest is purely human genes, after an intervention (inflammation-related genes).

About the pippeting issue. As you know is never discarded, but the internal controls or reference genes (as I am adressing them) goal is to normalize results and remove potential quantification errors. However, as you said the RNA concentration (A260) and purity were calculated spectrophotometrically from peripheral blood samples origin

As far as I know, the terms correlation and covariance related to PCA, are meant to discrimanate the standardization process (in correlation matrix) and not standardized matrix (scale function in R). Not sure if I understood this is what you imply?

My guess is nowadays everything is RNAseq and microarrays, whose mechanisms, and therefore, statistical preprocesses are different from RT-qPCR. But there are papers of PCA working wiht data coming from RT-qPCR. Shouldn't it work though?

**Update2**

Covariance and correlation: *"Covariance analysis of PCA: PCA is dimensionality reduction based on maximising multi-dimensional variance so that would be unusual, but correlation matrix of PCA is understandable. Normally covariance analysis is performed directly on standardised data."* Let me insist on this topic a little bit more.I am not sure if I understood correctly what you stated, but the information I've reviewed seems to say otherwise. As far as I was able to find, https://aedin.github.io/PCAworkshop/articles/b\_PCA.html : "*PCA was computed as a singular value decomposition (SVD) of a column centered, scaled matrix. This was PCA of the correlation matrix. If the matrix is centered but not scaled, it is PCA of the covariance matrix.*" I am not sure if I understood correctly what you stated, but the information I've reviewed seems to say otherwise

"*The correlation matrix is the standardized version of the covariance matrix. Analysing the correlation matrix is a useful default method because it takes the standardized form of the matrix; therefore, if variables have been measured using different scales this will not affect the analysis. Often you will want to analyse variables that use different measurement scales. Analysing the correlation matrix ensures that differences in measurement scales are accounted for. In addition, even variables measured using the same scale can have very different variances and this too creates problems for principal component analysis. Using the correlation matrix eliminates this problem also. There are statistical reasons for preferring to analyse the covariance matrix (the reason being that correlation coefficients are insensitive to variations in the dispersion of data whereas covariance is and so produces better-defined factor structures (Tinsley & Tinsley, 1987)) and generally the results will differ from analysis on the correlation matrix. However, the covariance matrix should be analysed only when your variables are commensurable.*" Conclusion: shouldn't it be the covariance matrix the one without standardized values, because variables share same scale?

Missing values: this is an issue, at least with several functions in R (prcomp(), princomp())which don't work with NAs in the matrix.

About this statement: "The key question, which is standardisation and I would say it is essential for RT-qPCR derived data at a minimum to at least investigate consistency of signal against a control. That minimum is what the second equation describes (why it works is complicated)."

I don't really understand what you mean standardize at a minimum. The endogenous control genes we are using are individualized for every sample. Let say we are investigating the fold-change in the individual **i** and the variable **k**. We are using endogenous controls in the **i** individual to "standardize" values. The

# Answer
> 0 votes
This is as good an answer I can give based the information. The issue here is that the experiment is derived via RT-qPCR and thats not the same as RNAseq type approach that the normal basis for questions on this site. Generally, we sort of assume RNAseq for gene expression where all the gene expression within a sample are automatically standardised against each other. Thus PCA can be performed directly on that data type. I will summarise this at the end because it's the key point I'm trying to make.

In your experiment there are 51 separate reactions per patient (or is 51 the total time series e.g. for one gene??) and it's not clear whether this constitutes a time series (but there must be some time-series aspect to the study: I assume the 'baseline' is at least one point in a time series). That is a lot of opportunity for variation within **one** patient sample, even if it is the same time point. Maybe the tech doing the RT-qPCR makes a pipetting error for 1 in 10 samples - the results could collapse. The other issue is whether the RT-qPCR is internally standardised, i.e. are you assessing two amplicons per PCR and the same control amplicon is used for each different gene under investigation? This traps the 1 in 10 hypothetical error rate because that is the first point of standardisation. The question is whether a further standardisation is required beyond that, i.e. per patient - thats tricky but basically its difficult to avoid because of variation between different samples from the same patient (again you must be doing some sort of time series). So for example the RNA extraction could have different amounts of blood etc ... However, it needs some thought because here you are standardisation a standardised result against a backdrop of an unstated hypothesis (comments).

**Summary** The single point I want to make is that large scale RT-qPCR is quite different from the standard RNA-seq approach on this site from a statistical perspective, i.e. it appear to be 51 reactions per sample whereas normally the approach is 1 reaction per sample which generates all the expression data across thousands of genes, thus PCA is a fantastic analysis for this data-type. This difference creates a layer of statistical complexity and I what I would do is carefully consider my hypothesis and formulate an analytical plan in accordance with the aims.

The question is whether a simple standardisation could be performed and then shift straight into PCA? Yes. Is that sufficient? I dunno. I think you might be avoiding time-series analysis by standardisation and then subtraction of the baseline which is then placed into PCA. At a minimum I would perform two PCA analyses one for baseline and one for post-treatment.

Generally, a complex standardisation (like I'm inferring) is not cool if the next stage of the analysis is PCA.

**Final thoughts** Although PCA is considered the primary analysis here. I think it takes serious thinking given all the available information to formulate an analytical plan from my experience in this type of data, but as a rough guide (given I dunno the experimental aims) a complex standardisation and then I personally would consider covariance thereafter for a time series rather than PCA as general approach here. *However*, two PCAs before and after using a simple standardisation and compare presence/absence, that could work and very clear results would be accepted.

---

**Note** Just to note that "correlation matrix based PCA" is basically PCA, it's just trying to sort out the groups. However "covariance" is a different analysis to PCA and statistically and the aims are very different. Also just to mention, but I assume you know this, "correlation" is different from "correlation matrix based PCA".

**Missing values** Missing values for PCA isn't an issue - you just miss them out. Missing values in a time series can have lots of different solutions, the key is justification.

---

**Update**

Covariance analysis of PCA: PCA is dimensionality reduction based on maximising multi-dimensional variance so that would be unusual, but correlation matrix of PCA is understandable. Normally covariance analysis is performed directly on standardised data.

Overall, I am now clearer one part of the calculation: this is essentially an *ad hoc* statistic to assess fold-change of the raw *Ct* via a single /minimal PCAs following one or more standardisations.

What you are measuring is the optimised multi-dimensional variance of the fold-change of raw *Ct* after at least one or else numerous standardisations. So this is not the actual fold-change of *Ct*, but variance of the fold-change derived via Eigen values.

If you are clear about that and what this represents and the reader is clear about this then it would be up to the reviewers.

The key question, which is standardisation and I would say it is essential for RT-qPCR derived data at a minimum to at least investigate consistency of signal against a control. That *minimum* is what the second equation describes (why it works is complicated).

---
Tags: r, gene-expression, pca
---


# Question
Title: How many percents explained variance by the first 50 principal components of PCA analysis is suitable for downstream `KNN` and `umap`?
In single-cell RNA sequencing analysis, I am wondering how many percents explained variance by the first 50 principal components of PCA analysis is suitable for downstream `KNN` and `umap`?

If there are lots of samples (\>40) and cells (\>100,000) to analyze together, there are a lot of `highly_variable_genes` (\>4000) can be found

```
import scanpy as sc
sc.pp.highly_variable_genes(adata, flavor='cell_ranger', n_top_genes=None)

```

The explained variance by the first 50 principal components of PCA analysis `sc.pp.pca(adata, n_comps=50, use_highly_variable=True, svd_solver='arpack')` is very low (\<12%), increase the number of principal components does not help a lot. I assume that the low explained variance of PCs will make `KNN` and `umap` difficult to separate the cells well. Am I correct?

```
sc.pp.neighbors(adata, n_neighbors=15)
sc.tl.umap(adata, min_dist=0.3, spread=2.0)

```

I can reduce the number of `highly_variable_genes` to get more explained variance by the first 50 principal components, but `KNN` and `umap` still can not separate the cells well.

I checked the first 5 PCs when I reduce the number of highly\_variable\_genes, and it doesn't look very high:

```
>>> np.cumsum(adata.uns['pca']['variance_ratio'])

array([0.08419679, 0.12635055, 0.13610176, 0.14490107, 0.15326007])

```

Any suggestions would be greatly appreciated.

# Answer
> 1 votes
**Edit** To answer the question, you're looking for 80% variance .... from the comments the actual command is `np.cumsum` which changes things a lot.

Previous response ... Firstly it's still good work. It wouldn't be sensible to proceed with `KNN` if 12% explained 50 PCs, now there are 15% in 5 PCs.

You wrote,

> 0.08419679, 0.12635055, 0.13610176, 0.14490107, 0.15326007

PC1 ~ 0.085, PC2 ~ 0.068

I think that you've cross checked the analysis and that is a good sign of good analysis.

The only easy possibility is to completely tighten down on all "variable genes" parameter. The total variance for the analysis think 50% total variance rather than an idealistic 80%.

The emergency measures ...

* Simply use say 10 samples
* Reducing the total number of genes, e.g. removing more variable genes looks good, anything else
* Check the input internally to make sure there is no formatting bug
* Check via Seraut

The thing I would cross check is reducing the sample size a lot ... is there a stray sample(s) here that is throwing everything out? That would be thing I'd really look into. What happens to the total variance when say 20 samples, or 10 samples are used?

Just to finalise, are there any natural population divisions in the sample that could be used? If this results in a huge jump in variance that would be an answer. The only thing that is unusual is that there are a lot of samples for this type of study (from my perspective).

---
Tags: single-cell
---


# Question
Title: Why are conversion output SMARTS strings are not the same for different file types of the same structure?
I've just used *Open Babel* to convert *.sdf* and *.mol2* files to *SMARTS*, and I'm a bit surprised the output can be different for the same structure. It's confusing me that one structure, one piece of information can be different depending on the file format in which it was saved.

I converted a ligand *1N7* in two available formats (*.sdf* and *.mol2*) of the *7RDX* structure from *RSCB PDB* for chain *N*, and also for chain *U*. After the conversion, the *SMARTS* strings for the ligand in both chains are not the same for both input files types (i.e. the difference between the same ligand, from the same chain, just in two different formats available to download from RSCB PDB):

# example no. 1:

\[additional `@H` added in places, some square brackets added\]

```
$ cat ./7rdx_N_1N7_from_mol2.smi

C1[C@@]2(C)[C@@H](C[C@@H](C1)O)[C][C]([C]1[C@H]3[C@]([C@H](C[C]21)O)([C@H](CC3)[C]([C])[C][C])C)O   1N7

$ cat ./7rdx_N_1N7_from_sdf.smi

C1[C@@]2(C)[C@@H](C[C@@H](C1)O)C[C@H]([C@H]1[C@H]3[C@]([C@H](C[C@H]21)O)([C@H](CC3)[C@H](C)CC)C)O   1N7

```

# example no. 2:

\[some square brackets added\]

```
$ cat ./7rdx_U_1N7_from_mol2.smi

C1[C@@]2(C)[C@@H](C[C@@H](C1)O)C[C@H]([C@H]1[C@H]3[C@]([C@H](C[C@H]21)O)([C@H](CC3)[C@H](C)CCC(=O)NC[C][C][N+]([C])([C])[C])C)O 1N7

$ cat ./7rdx_U_1N7_from_sdf.smi

C1[C@@]2(C)[C@@H](C[C@@H](C1)O)C[C@H]([C@H]1[C@H]3[C@]([C@H](C[C@H]21)O)([C@H](CC3)[C@H](C)CCC(=O)NCCC[N+](C)(C)C)C)O   1N7

```

In both examples, the structures look quite similar to each other, but the difference is a little worrisome, because there is no chance to check every output by hand.

Is it a normal behaviour during the conversion?
Is it related to a different notation in the input files?
Which input format is better?

# Answer
to me your 2 strings seem different, no N in first molecule.

Here pictures, I used a PDB version of the original mmCIF file

but dindt use the RCSB PDB available one, imported .cif into

Biopython and parsed the resulting Struct object \[so there could be mistakes\].

so chain N --\> A and chain U --\> E

as per in my comment:

> Mmmhhh be sure that each of the 3 modeled ligands correspond to the full one, could be that the density is poor and pdb shows only what is really ((??) there and not what they soaked the crystal with, actually 1N7 is CHAPSO I believe a detergent

... I could be wrong though, this mmCIF format is driving me crazy
> 1 votes

---
Tags: pdb
---


# Question
Title: Genomic inflation factor distribution (lambda)
Considering that the p-value of a chi statistic(df=1) is uniform, the expected genomic inflation factor is λ=1. I show its calculation in R with the next chunk of code:

```
chi = qchisq(df$pval,df=1,lower.tail=FALSE)
factor = median(chi)/qchisq(p=0.5,df=1)

```

My question comes here. Imagine that I compute λ on several (e.g. \>30 independent GWAS) studies. My null hypothesis is that all of them lack true-polygenic signals (i.e. uniform p-value distribution). What test can I do to reject H0? What is the distribution that I can assume for the inflation factors?

# Answer
The distribution as I suspected seems very much like a normal distribution. However if someone knows to explain the reason (limit central theorem?) I will give select their answer.

```
    # Do simulation
    factor_calc = function(pvals){
      chi = qchisq(pvals,df=1,lower.tail=FALSE)
      return(median(chi)/qchisq(p=0.5,df=1))}
    pvals_studies = sapply(1:10000, function(i){runif(10000)})
    factor_distribution = apply(pvals_studies, 2, function(pvals){factor_calc(pvals)})

    # Plot resulting distribution
    hist(factor_distribution)

```
> 2 votes

---
Tags: statistics, gwas
---


title: Display row and column names in R

text:
# Question
Title: Display row and column names in R
I have plotted a heat-map with 235 rows and 570 columns which is below: how do we make clear names of row and columns on this graph using R?

# Answer
You don't, simple as that, it is too many in both columns and rows to show them all. You could shrink size, but then it becomes unreadable. You can highlight individual genes (rows) that are important to support the message you want to send with that experiment, with packages such as ComplexHeatmap, see the Mark annotation section.

Alternatively, you can use something like hierarchical clustering to form groups, then make functional annotations (such as REACTOME terms) and then annotate eaqch cluster with selected terms, again to support the message you want to send. Go through the ComplexHeatmap manual, lots of examples with code there to get inspiration from.

As for the columns, you could use block annotations to color-highlight columns bhelonging to the same experimental group, or split the heatmap by these groups to put some visual separation/guidance.
> 4 votes

# Answer
*\[note: This is a generic answer; the post was not specific enough for me to be able to include any more useful detail\]*

A heatmap of 235 rows and 570 columns is within the realm of what should be explorable / visible with a bit of tweaking. I have needed to do this with large heatmap plots - for example, single cell data with 500 genes and 2000 cells, where I care about the genes but not the cell IDs. My usual approach is to make the dimensions of the output image larger, and also reduce the text size.

The method of doing that depends a lot on how the plot was generated initially. For example, if using a base R package:

```
pdf("output_file_name.pdf", width=57, height=23.5); # about 10 squares per inch
## <plot generation function>
dev.off()

```

If using ggplot:

```
library(tidyverse);
input_data %>%
  ggplot() +
  aes(...) +
  # <plot generation function>
ggsave("output_file_name.pdf", width=57, height=23.5);

```

Changing the function used to display the heatmap can sometimes help as well (e.g. `pheatmap` instead of `heatmap`).

Within the functions, text size can often be reduced to help with this as well (e.g. `cex = 0.2`). There is a limit to how much this can be done, as most PDF viewing programs have a zoom limit - I think Adobe Reader's limit is 6400% (i.e. 64 times magnification) - as well as maximum page dimension limits.
> 2 votes

---
Tags: rna-seq, ggplot2
---


# Question
Title: How to apply a mathematical function to a matrix in R
I have a matrix that looks something like this:

| chr | pos | het |
| --- | --- | --- |
| chrI | 27 | 2 |
| chrI | 55 | 0 |
| chrI | 27 | 0 |
| chrI | 55 | 1 |

It's essentially a VCF file that I have converted into a TSV. I want to be able to take all the values in the `het` column and apply this function: `het_df['het'] = het_df['het'].transform(lambda x: -1*(x-1))`

That line is from a python script and it works really well in transforming that column. However, as I am making plots from this data I want to be able to do this in R. When I went to convert it to an R script I started with this function:

`convert_df <- lapply(het_df[,3], FUN=function(x) formula=(-1*(x-1)))`

As I saw it in this issue: https://stackoverflow.com/questions/7009333/how-to-use-lapply-with-a-formula

However, I keep getting error messages such as this one:

```
Error in x - 1 : non-numeric argument to binary operator

```

Which is odd because when I run `het_df[,3]` it gives me the correct output:

```
   [1] "2" "0" "1" "1" "1" "1" "1"

```

I think I'm just a little lost and I don't know what to do to run this function...everything worked fine in my python script (I was able to generate plots with matplotlib) but I like manipulating plots more in R than python so I wanted to be able to convert the script. If anyone is able to see what I'm clearly doing wrong please let me know! Thank you!

Here is the python code I've been using in case this helps (I'm transforming the last column so I always get a 1 or a 0 for heterozygosity/homosygosity and then using a rolling window to obtain the variability in the data series and plotting them):

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

# read in file
het_df = pd.read_csv("het_file.tsv", sep="\t", header=None).set_axis(['chr', 'pos', 'het'], axis=1, copy=False)

# getting rid of homozygous alt allele
# always get 1 for heterozygosity
het_df['het'] = het_df['het'].transform(lambda x: -1*(x-1))

windows = [25, 50, 100, 200] # 100 bp, get rid of a lot of noise
steps = [5]

for window in windows:
    for step in steps:
        het_ratio = {}
        for chr_het_df in het_df.groupby('chr'):
            # rolling over the window
            # center = have the count in the center of the window
            # 50bp window count in the front or back of window or count in the center
            # ratio of heterozygosity
            het_ratio[chr_het_df[0]] = chr_het_df[1].set_index('pos').rolling(window=window, step=step, center=True).mean()
        fig, ax = plt.subplots(len(het_ratio), 1, sharex=True)
        for index in range(len(het_ratio)):
            chromosome = list(het_ratio.keys())[index]
            sns.scatterplot(het_ratio[chromosome], ax=ax[index], s=4)
            ax[index].set_ylim(bottom=0, top=1)
            ax[index].legend([], [], frameon=False)
            ax[index].set(xlabel="Chromosome position", ylabel=chromosome)
            ax[index].set_yticks([0, 0.5, 1])
            ax[index].set_yticklabels([0, 0.5, 1], size = 6)
        fig.subplots_adjust(top=0.98, bottom=0.05, left=0.05, right=0.98)
        fig.set_size_inches(20, 15)
        plt.savefig("het_plot" + str(window) + "_step" + str(step) + ".png", dpi=150)
```
</code></pre>

# Answer
If the numbers are in quotes it means that it has read them as text. A numeric vector will be:

\[1\] 2 0 1 1 1 1 1

To replace that column with the new value you can use this:

```
het_df[,3] <- sapply(as.numeric(het_df[,3]), FUN=function(x){-1*(x-1)})

```

I used `het_df[, 3] <-`  to replace the content of the column. I changed `lapply` to `sapply`, to get a vector and not a list. I forced the input to be numeric via `as.numeric`.
The syntax of the functions was wrong, to create a function in R you need `function(x){#code}`. The linked question is for a function that has an argument called formula, which is not the case here, because you want to define your own (anonymous) function.
> 1 votes

# Answer
Doing mathematical operations on a matrix will only work if *all* the matrix elements are numerical. In your case, you have mixed character and numerical variables, and it looks like a data frame or tibble would be more appropriate (calling it `het_df` is a strong indicator that this is actually what you want to do, as well as the `pandas` structures in the python code). You can load a TSV file into a tibble in R using the `read_tsv` function from `tidyverse`:

```
library(tidyverse);
het_df <- read_tsv("input_file_name.tsv");

```

After this, numerical mathematical operations can be carried out column-by-column, because the columns are treated as separate structures internally:

```
het_df$het <- -1 * (het_df$het - 1); # approach 1 [base R]
het_df <- mutate(het_df, het = -1 * (het - 1); # approach 2 [tidyverse]

```
> 0 votes

---
Tags: r, python
---


# Question
Title: Merge the dataframes with python iterator
I have 3 experiment files (`File1.tsv, File2.tsv, File3.tsv`) from three experiments and 1 library file (`library_file.tsv`). I want to merge the files by creating a python iterator because the number of experiment files will increase with every additional experiment and the library file will be always the same. All the experiment files are in the same directory.

`File 1`

```
id  gene    Exp_1
ID000000001 EGFR    0.127
ID000000002 ADAD1   0.421
ID000000003 ADAD1   0.301

```

`File2`

```
id  gene    Exp2
ID000000001 EGFR    0.303
ID000000002 ADAD1   0.236
ID000000003 ADAD1   0.071

```

`File3`

```
id  gene    Exp_3
ID000000001 EGFR    0.847
ID000000002 ADAD1   0.992
ID000000003 ADAD1   0.912

```

`library_file`

```
id  value1  value2
ID000000001 Single  Guides_Single
ID000000002 Single  Guides_Single
ID000000003 Single  Guides_Single

```

`Expected output`

```
id  gene    Exp_1   Exp2    Exp_3   value1  value2
ID000000001 EGFR    0.127   0.303   0.847   Single  Guides_Single
ID000000002 ADAD1   0.421   0.236   0.992   Single  Guides_Single
ID000000003 ADAD1   0.301   0.071   0.912   Single  Guides_Single

```

# Answer
> 2 votes
Its a simple join operation.

```
import pandas as pd

df1 = pd.DataFrame({'id':['ID000000001','ID000000002','ID000000003'], 'gene1':['EGFR', 'ADAD1','ADAD1'],'score1':[0.127, 0.421, 0.301]}).set_index('id')
df2 = pd.DataFrame({'id':['ID000000001','ID000000002','ID000000003'], 'gene2':['EGFR', 'ADAD1','ADAD1'], 'score2':[0.303, 0.236, 0.071]}).set_index('id')
df3 = pd.DataFrame({'id':['ID000000001','ID000000002','ID000000003'], 'gene3':['EGFR', 'ADAD1','ADAD1'], 'score3':[0.847,0.992,0.912]}).set_index('id')
df4 = pd.DataFrame({'id':['ID000000001','ID000000002','ID000000003'], 'v1':['Single', 'Single', 'Single'], 'v2':['Guide', 'Guide', 'Guide']}).set_index('id')

df = df1.join([df2,df3,df4])
df.drop(['gene2', 'gene3'],axis=1,inplace=True)
print (df)

```

Output

```
             gene1  score1  score2  score3      v1     v2
id
ID000000001   EGFR   0.127   0.303   0.847  Single  Guide
ID000000002  ADAD1   0.421   0.236   0.992  Single  Guide
ID000000003  ADAD1   0.301   0.071   0.912  Single  Guide

```

'Iterator' just `pickle` via

```
path='/home/desktop/df.pkl'
df.to_pickle(path)
# Reload
df = pd.read_pickle(path)
# Add another data set
df = df.join(newdf)
df.drop(['gene4'],axis=1,inplace=True)

```

... then pickle it again

---
Tags: python, pandas, merge
---


title: Sanger sequencing annotation error

text:
# Question
Title: Sanger sequencing annotation error
I am a student in a Cancer lab. Working with sanger is new to me. While analyzing a report we found an insertion that has not been reported in any databases so far, we were working on checking if the insertion is pathogenic or not but when we reanalyzed the sanger annotation through variant reporter of the insertion it was giving it one position ahead and we analyzed that with the VEP as well that is were the confusion started.

The sequencer showed annotation of an insertion not between the bases but a step/position next to it:

**c.5532-5533insGAA** is the insertion position came after position 32 and before 33 but the sequencer annotates it like **c.5533-5534insGAA**.

We are working on human genome hg38.

Should I consider it a machine error or should I take that position for the result? Or is there a nomenclature error I am making here?

And the annotation is in this form \[c.33-34insAUU\]+\[=\] and I don't know what the +\[=\] represents here.

# Answer
Just to give my thoughts:

* The missing information is the actual sequence

1. Is this a conflict in the change in position of the insert?
2. Could it simply be a consequence of down stream numbering, either a logical error (coding stuff), or lab artefact.

My bet is its point 2, but really we don't know except for the OP because we don't know the sequence.

**Sources of error, Point 2**

* Sanger sequencing was traditionally known for "frame-shift" errors if there were mono-mers, e.g. AAAA -\> AAAAA, it can mess up over the number of bases in a mono-mer.
* PCR amplification can also cause slippage against monomers, same thing as the above point - but it's real.
* 1 vs 0 ... this is because in a Pythonic list or any C/Perl-based array numbering starts at 0 (not 1) - hence @terdon's comments.

If upstream mono-mers are present that might be the source of the confusion here and given the amount of sequence data thats very likely.

What I'd do:

* Just Blast the new sequence and see what emerges: if it is genuinely novel - they'll be zero 100% identity hits.
* VCFs ... I think this is VarSome's stuff (forget).

If there are Blast hits it will identify whether there is any upstream lab artefacts.

**Conclusion** If this was any other situation in human genetics, or even mammalian genetics I'd be excited. However, in this case - I think its a routine result. Firstly, the chance of finding a new "driver" mutation are slim. Secondly, I've looked at this chromosome in a cancer patient (only once) and it was nuts (and thats an under statement) - I've a fairly good memory for numerical patterns, and I'd never seen a pathogenic species behaving like that. I think this insert is just a consequence of the entire mutational chaos that are cancer cells.

Just to make it clear: I don't do cancer but if you want an authoritative assessment and think you've something major I'd contact (or your PI) the Broad Institute via Geraldine A. Van der Auwera who also has an account on this site via @Geraldine\_VdAuwera. Sanger sequencing is not what they do, but they'd know at the drop of a hat what the significance and if they don't respond - thats your answer ;-). Alternatively, VarSome would know, albeit they're commercial enterprise offering freeware.
> 0 votes

---
Tags: ngs, sequence-annotation, assembly, mutations, sequence
---


# Question
Title: How to change the position of groups colors in DimPlot of Seurat
I have a Seurat object and plotted the Dimplot for UMAP visualization for 2 variables, as shown in the image below.

Now, the problem is that I want the group by variables such as `Non-responder and Responder` and `anti-CLTA4, anti-CLTA4+PD1, anti-PD1` on the top of the UMAP plot and not on the right side. My desired output would look like the following.

The timely help is highly appreciated, many thanks.

# Answer
DimPlot returns a ggplot object, so ggplot functions can be applied to it. You can change the legend position like this:

```
DimPlot(data) + theme(legend.position = "top")
# or
DimPlot(data) + theme(legend.position = c(.1, .9))

```
> 3 votes

---
Tags: scrnaseq, seurat, ggplot2
---


# Question
Title: Interpreting short indel calls in 1000 Genomes Project VCFs
Consider the following short indel polymorphism rs59679400 on chr7.

```
$ rsidx search dbSNP/dbSNP_GRCh38.{vcf.gz,rsidx} rs59679400 | cut -f 1-7
[rsidx] running version 0.3
NC_000007.14    12510038        rs59679400      GAAAA   G,GA,GAA,GAAA,GAAAAA    .       .

```

In the Phase III (2015) 1000 Genomes VCFs (GRCh37), I see the following record for the variant.

```
$ rsidx search 1000GenomesPhaseIII/ALL.chr7.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.{vcf.gz,rsidx} rs59679400 | cut -f 1-7
[rsidx] running version 0.3
7       12549664        rs59679400      GA      G       100     PASS

```

That's easy enough to interpret, although it looks like only two of the several possible alleles at that locus are represented.

When I look at the VCFs from the latest NYGC update to the 1000 Genomes data, I see the following. (*I'll note that it's really inconvenient that RS numbers are no longer included for any variants.*)

```
$ tabix 1000GenomesNYGC/1kGP_high_coverage_Illumina.chr7.filtered.SNV_INDEL_SV_phased_panel.vcf.gz chr7:12510038-12510038 | cut -f 1-7
chr7    12474780        HGSV_106040     G       <DUP>   .       .
chr7    12510038        7:12510038:GA:G GA      G       .       .
chr7    12510038        7:12510039:GAA:G        GAA     G       .       .

```

Ignoring the first record, which is a larger structural variant, we see two records at the position in question. **My initial question is why they've used two records to encode variation for a single indel variant?** As shown in the dbSNP record in the first snippet, VCF allows one to encode an arbitrary number of alternate alleles in a single record.

**My second—and more important—question is how one is to interpret these records and resolve the sequence for any particular individual genotyped at the locus?** Selecting an arbitrary sample `HG00103`, I get the following allele calls for the two records? What am I to make of this?

```
>>> from pysam import VariantFile
>>> vcf = VariantFile("1000GenomesNYGC/1kGP_high_coverage_Illumina.chr7.filtered.SNV_INDEL_SV_phased_panel.vcf.gz")
>>> records = list(vcf.fetch("chr7", 12510037, 12510038))
>>> len(records)
3
>>> records[1].samples["HG00103"].alleles
('G', 'G')
>>> records[2].samples["HG00103"].alleles
('GAA', 'GAA')

```

# Answer
> 2 votes
this is Marta from the NYGC. Thanks for reaching out with your questions.

The full SNV/INDEL callset that we released (listed as "SNV/INDEL VCFs" in the key resource table in our paper) includes multiallelic representation of all multiallelic sites (i.e. multiple ALT alleles are listed together in column 5 separated by commas) and the record you showed as an example above is listed there in the following way:

chr7 12510038 . GAAA GAA,GA,G

The VCF file you were looking at is the phased panel that we generated and released alongside the full unphased callset. The phasing software tools that we used in the paper (SHAPEIT2, SHAPEIT4, Eagle2) as well as most other phasing tools do not support multiallelic sites and require that such sites are split into separate lines and represented as biallelic for the purposes of phasing. As you noted, this is not ideal and can lead to difficulties in interpretation (as others explained above) but it's the best we can do with the existing approaches. We mention this in the Methods section of the paper as well as in the README that's included with the phased panel on the IGSR ftp (all linked in the paper). Please let us know if you have any further questions related to the high-coverage 1kGP resources. Thanks!

# Answer
> 3 votes
Older GATK was using the so-called variant model where overlapping variants were combined during calling. More recent GATK is using the allelic model where one non-reference allele is evaluated at a time.

With the allelic model, a reference allele doesn't mean the sample possesses a reference allele; it means the sample doesn't have the alternate allele. For this particular example, records\[1\] means a homozygous -A deletion. records\[2\] means the sample doesn't have the -AA deletion. The two records are consistent with each other.

It is still possible to have conflictive records (for example, homozygous -A deletion and homozygous -AA deletion at the same location). At least this used to be a problem with GATK. Not sure if/how much this has been improved recently.

# Answer
> 1 votes
Often multi-allelic sites would be split into several lines (one per alternate allele) to make filtering easier for specific types of alternate alleles (e.g. filter out MNV vs SNV).

VCF-level, one thing you can do to identify the actual genotype over the locus easier is to join the split alleles into a single multi-allelic record. You can easily do this using bcftools:

```
bcftools view --samples HG00103 --regions chr7:12510038-12510038 1kGP_high_coverage_Illumina.chr7.filtered.SNV_INDEL_SV_phased_panel.vcf.gz | bcftools norm -m + | bcftools view -H

```

From here is clear that starting on position chr7:12510038, instead of "GAA", you have homozygous "GA" (1|1), therefore a hmozygous "A" deletion, which could be either the first or second "A" at positions 12510039 and 12510040.

---
Tags: vcf, 1000genomes
---


title: Compare my VCF to gnomAD variants

text:
# Question
Title: Compare my VCF to gnomAD variants
I have a VCF with small variant calls against HG38 and I would like to determine which of those calls are present in the gnomAD database.

Is there an existing tool that can do this? Should I be looking at variant annotation tools? or is this something that I would look at BCFtools for?

# Answer
> 3 votes
You can use vcfanno to annotate your vcf with the gnomad vcf file. There's an example config file for annotating with gnomad (along with other databases such as dbSNP) here.

# Answer
> 2 votes
I haven't checked how fast this is but you could try using Hail (disclaimer: I'm a maintainer). You'll need to install the GCS connector first to read from Google Cloud Storage:

```
curl https://broad.io/install-gcs-connector | python3

```

```
import hail as hl

mt = hl.import_vcf('your.vcf')
mt = hl.split_multi_hts(mt)  # use the same biallelic representation that gnomAD uses
mt.write('data.mt')  # store in Hail's native format for faster queries
mt = hl.read_matrix_table('data.mt')
sites = hl.read_table(
    'gs://gcp-public-data--gnomad/release/3.1.2/ht/genomes/gnomad.genomes.v3.1.2.sites.ht'
)
mt = mt.anti_join_rows(sites)

```

That'll filter down to just those variants (defined by chromosome, position, and alleles) which are not present in gnomAD v3.

# Answer
> 1 votes
If you have sufficient available disk space, you could download the files from the gnomad downloads page and use a bedtools intersect command to find the overlapping variants.

However, the recommended platforms for gnomAD are to use hail utilities developed by them and others, due to the constraints of working with ~1TB of data.

I am not an expert but potentially someone else knows more about using the site hail data.

---
Tags: vcf
---


# Question
Title: Simulating phenotype with the 1000 Genomes Project
I'm looking for a way to simulate phenotypes against a real SNP data source, such as the 1000 Genomes. It must be free for commercial purpose (Eg.: MIT license). Any recommendation? I'm trying to use the GCTA64, but I couldn't get it working. The documentation doesn't help much. At the end of the day, I want: 1 - Simulate Case/Control and/or quantitative phenotype 2 - Link it with a real SNP dataset (eg: 1000 Genomes) 3 - Conduct GWAS analysis using Plink and/or Hail.

Looking Hail's and Plink's tutorials for GWAS, I realised both use simulated phenotypes from real SNP data sources (1000 Genomes and HapMap, respectively), but how they created the datasets are beyond of the scope of the tutorials, thus not reported.

As mentioned before, I've tried `gcta64`, but no success. Here's what I've tried:

1 - Downloaded 1000 Genome sample from Plink page: Entire dataset as a single .tar.gz (1.12 GB) (A2 allele major, not ref, on chr3 before 15 Oct 2017)

2 - Tried to generate the simulate data by:

```
./gcta64  --bfile 1kGenomesP1/1kg_phase1_all/1kg_phase1_all --simu-qt  --simu-causal-loci causal.snplist  --simu-hsq 0.5 --simu-rep 3  --keep test
.indi.list --out 1kg_phase1_all

```

Error: `Error: --keep test.indi.list not found`.

Question: What should be the files:`causal.snplist`, `test.indi.list`, `teest`? I know the property name suggests what should be, but there's no example of file structure nowhere... Any practical example or tutorial?

Btw - I apologise in advance if this is a too trivial question. I'm quite new at it. Appreciate your patience and help :)

# Answer
I managed to get it done without the gcta64. I have written a Jupyter Notebook with the steps: https://github.com/bambrozio/bioinformatics/blob/master/utils/1k-genomes-phenotype-simulator.ipynb
> 1 votes

# Answer
You might try `ldscsim` which was developed by Nik Baya for this purpose.

https://hail.is/docs/0.2/experimental/ldscsim.html

```
import hail as hl

mt = hl.import_vcf('path/to/your/data.vcf')
mt.write('data.mt') # store in Hail's native format for speed
mt = hl.read_matrix_table('data.mt')
mt = hl.experimental.ldscsim(
    mt, mt.GT, h2=0.5
)

```

`ldscsim` can generate phenotypes under a number of different models.

For the GWAS tutorial, AFAIK, we just simulated betas for each site, then computed the phenotype from that. If you ask on https://discuss.hail.is and tag Tim Poterba he can probably give you the precise answer.

```
def logistic(x):
    return 1 / (1 + hl.exp(-x))

mt = mt.annotate_rows(beta = hl.rand_norm(0, 1))
mt = mt.annotate_cols(
    is_case = 0.5 < logistic(hl.agg.sum(beta * mt.GT.n_alt_alleles()))
)

```
> 1 votes

---
Tags: snp, gwas, 1000genomes
---


# Question
Title: If fastp output is not a good measure of FASTQ correctness, what is?
In the beginning of my pipeline, I just fed the paired reads (2 files) into fastp, with the default options, and assumed it would do a good job preparing the reads for the next step: alignment

But I now checked the files in FastQC and it doesn't show much difference from before the fastp processing: Some extra yellow and red signals appear. I think the end result is even worse.

fastqc doesn't analyse paired reads files together. So it might be wrong regarding the paired quality, taken as a whole.

Anything else I should use together with FastQC, to verify the output quality?

I'm aligning a bunch of poolseqs and then comparing them to find consistency in the change of snp frequencies

As requested here is one of the examples that made me ask this:

***Before fastp. only one problem noted by fastqc***

***After fastp. overall quality seems to have decreased***

The raw data is a pool seq from illumina.

# Answer
> 3 votes
On a regular basis, using FastQC is quite enough as to assess the quality of FastQ-formatted data. It gives you plenty of details as can be seen in the screenshots you shared.

As others have commented or suggested, there is nothing wrong with having more 'yellow' warnings after processing the FASTQ files using fastp. At lest not as long as what you observe in the data looks fine and the warning has some explanation.

The sequence length distribution warning is very often present post-processing since at the begging all your reads are of the same length but this is not the case any more after you trim the adapters and the low quality bases at the end of the read for example. You have to go to the plots specific to that section and make sure all looks fine. It should be the case.

The GC sequence content warning. It depends a lot what organism you are sequencing, which regions of the genomes, etc. Importantly, what is the expected GC content of those regions. So you can compare to these observed numbers. Once again the warning can be easily triggered after the processing, since for a start Illumina sequencing has known GC bias

Finally, I think actually your reads look better after the processing with fastp. You are showing us the per base sequence quality and before, there are many reads with bases on the yellow and even the red sections of the quality scale. But this is not the case any more after the processing !

---
Tags: fastq, reads, fastqc, read-correction
---


title: What is changed in this PDB file?

text:
# Question
Title: What is changed in this PDB file?
I used a tool called reduce to add Hydrogen bonds information to the PDB file.

I used the command:

```
$ ./reduce -NOFLIP 2GB1.pdb > 2GB1_with_H.pdb

```

The output file is supplied as follows:

I haven't understood the output.

What special change was made by the tool to the `2GB1.pdb` file?

# Answer
> 5 votes
So the 'right' answer here would be for you to use a `diff` tool of some kind (WinDiff, or just the diff command in a Unix-like environment). That would show you the exact added lines. (Although, trying this myself does not produce very clear results...)

However, I know from running `reduce` on PDB files that it's added a few things:

* The "USER MOD" lines, there's one on the first line and a large block just before the ATOM records
* The hydrogen atoms, of course, needed for correct H-bond calculations

As some commenters have pointed out, this should really be in the paper or other documentation for the tool.

---
Tags: proteins, pdb
---


# Question
Title: Extract lines from a file based on a range of coordinates listed on a second file
I would be highly appreciated to have a script able to perform the following operation:

Only extract scaffold\_1 lines from file1.vcf in which POS number (column 2) falls within column 2 and column 3 numbers of scaffold\_1 lines from listfile.txt (inclusive). Similarly, only extract scaffold\_2 lines from file1.vcf in which POS number falls within column 2 and column 3 numbers of scaffold\_2 lines from listfile.txt. In reality the vcf I have continues until scaffold\_12 and has 17,000 lines.

file1.vcf

```
#CHROM      POS ID  REF ALT QUAL    FILTER  INFO    FORMAT          COUNTS
scaffold_1  160 .   C   A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_1  210 .   AC  A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_1  340 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_1  470 .   T   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_1  580 .   A   AC  40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_1  630 .   CA  C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  40  .   C   A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  100 .   AC  A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  370 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  600 .   T   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  800 .   A   AC  40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2      850     .       C       CTA     40      PASS    .       GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  1120    .   CA  C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40

```

listfile.txt

```
scaffold_1  50  100
scaffold_1  200 350
scaffold_2  350 450
scaffold_2  750 850

```

The desired output file is the following:

desired\_file.vcf

```
#CHROM      POS ID  REF ALT QUAL    FILTER  INFO    FORMAT          COUNTS
scaffold_1  210 .   AC  A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_1  340 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  370 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  800 .   A   AC  40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2      850     .       C       CTA     40      PASS    .       GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40

```

To do this job I found the below awk code from another similar post but they don't differentiate scaffold\_1 from scaffold\_2, and I don't how to add this extra variable properly. Thank you.

```
awk 'NR==FNR{rng++;start[rng]=$2;end[rng]=$3;next} {for (i=1;i<=rng;i++) if (($2>=start[i])&&($2<=end[i])) {print; next}}' listfile.txt file1.vcf > output.vcf

```

```
<listfile.txt awk '{ print "$2 >= " $2 " && $2 <= " $3 }' | awk -f - file1.vcf > output.vcf

```

I tried to modify the last awk code like this, but it didn't work:

```
<listfile.txt awk '{ print "$1 == " $1 " && $2 >= " $2 " && $2 <= " $3 }' | awk -f - file1.vcf > output.vcf

```

# Answer
Using GNU awk for arrays of arrays:

```
$ cat tst.awk
NR==FNR {
    ranges[$1][$2 FS $3]
    next
}
FNR == 1 {
    print
    next
}
$1 in ranges {
    for ( range in ranges[$1] ) {
        split(range,r)
        if ( (r[1] <= $2) && ($2 <= r[2]) ) {
            print
            next
        }
    }
}

```

```
$ awk -f tst.awk listfile.txt file1.vcf
#CHROM      POS ID  REF ALT QUAL    FILTER  INFO    FORMAT          COUNTS
scaffold_1  210 .   AC  A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_1  340 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  370 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  800 .   A   AC  40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2      850     .       C       CTA     40      PASS    .       GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40

```

Using any awk it'd be:

```
$ cat tst.awk
NR==FNR {
    ranges[$1] = ranges[$1] $2 FS $3 RS
    next
}
FNR == 1 {
    print
    next
}
$1 in ranges {
    n = split(ranges[$1],list,RS)
    for ( i=1; i<n; i++ ) {
        range = list[i]
        split(range,r)
        if ( (r[1] <= $2) && ($2 <= r[2]) ) {
            print
            next
        }
    }
}

```
> 4 votes

# Answer
If you have a VCF file with an intact header, you can use for this. All you need to do first is to bgzip and index your VCF:

```
$ bgzip file.vcf
$ bcftools index file.vcf.gz

```

Then, assuming your list of regions are actually 1-based tab-delimited, use `bcftools view` and the `--regions-file` option to restrict the search to the regions listed. Use also the `--no-header` option to suppress the header in the output:

```
bcftools view -HR regions.txt file.vcf.gz

```

From the common options documentation:

> -R, --regions-file FILE
>
> Regions can be specified either on command line or in a VCF, BED, or tab-delimited file (the default). The columns of the tab-delimited file can contain either positions (two-column format) or intervals (three-column format): CHROM, POS, and, optionally, END, where positions are 1-based and inclusive. The columns of the tab-delimited BED file are also CHROM, POS and END (trailing columns are ignored), but coordinates are 0-based, half-open. To indicate that a file be treated as BED rather than the 1-based tab-delimited file, the file must have the ".bed" or ".bed.gz" suffix (case-insensitive). Uncompressed files are stored in memory, while bgzip-compressed and tabix-indexed region files are streamed. Note that sequence names must match exactly, "chr20" is not the same as "20". Also note that chromosome ordering in FILE will be respected, the VCF will be processed in the order in which chromosomes first appear in FILE. However, within chromosomes, the VCF will always be processed in ascending genomic coordinate order no matter what order they appear in FILE. Note that overlapping regions in FILE can result in duplicated out of order positions in the output. This option requires indexed VCF/BCF files. Note that -R cannot be used in combination with -r.
> 7 votes

---
Tags: vcf, awk
---


title: Problem with ClinCNV script

text:
# Question
Title: Problem with ClinCNV script
I have some problems with ClinCNV:

My input: `Rscript /root/ClinCNV/clinCNV.R --normal /root/samples/final.cov --out /root/samples/ -bed /root/samples/fileGC.bed --folderWithScript /root/ClinCNV/`

And there is an output: `/root/ClinCNV/clinCNV.R: error: Error in getopt(spec = spec, opt = args) : short flag "b" is invalid`

I've already ask a question about ClinCNV and the main idea was to run script with absolute pathways.

# Answer
> 3 votes
You are using `-bed`. Multi-letter flags are usually set using `--` and not `-`. I don't know this particular program but since all other long options use `--`, and the only option with a `b` is your `-bed`, I suspect you just want `--bed` instead:

```
Rscript /root/ClinCNV/clinCNV.R --normal /root/samples/final.cov \
   --out /root/samples/ \
   --bed /root/samples/fileGC.bed \
   --folderWithScript /root/ClinCNV/

```

---

This is confirmed by the tool's documentation, actually:

> --bed - .bed file with coordinates of on-target regions, GC-annotated. Gene anotations in the last column can be used in output.

So yes, you need `--bed` and not `-bed`. Because you used `-`, the tool reads that as a single character flag (`-b`) with the value `ed` and doesn't know what to do with it.

---

On an unrelated note, you seem to be running these things as root which is a very, very bad idea. One bug in the program or one misconfiguration and you can mess up your entire system. I strongly urge you to never run commands as root unless absolutely necessary.

---
Tags: bed
---


title: Lifting snp rsid build 37 to 38

text:
# Question
Title: Lifting snp rsid build 37 to 38
I’m looking to lift over a dataset of GWAS summary stats from build 37 to 38, including converting rsid. Downstream it is to be used for co-localisation with another dataset of build 38.

Would anyone be able to help with code for doing this please?

# Answer
There are multiple tools to use for this. Might I suggest the following reads: Introduction to liftover tools, where they also refer to this usefull post BioStars post. Maybe also consider reading this paper for a word of caution: Paper.

Maybe you can also clarify you question more, if these links do not help you with what you need.
> 1 votes

# Answer
I have recently worked on a couple of tools to do exactly this. They are +munge to convert summary statistics into files following the GWAS-VCF specification and +liftover to convert a VCF from one genome build to another. This will automatically take care of REF⇆ALT switches by flipping the effect size and it will correctly handle indels, if there are any, in your summary statistics file
> 0 votes

---
Tags: gwas, identifiers, liftover
---


title: What do you use the CLANS program for?

text:
# Question
Title: What do you use the CLANS program for?
I am reading this paper called Birth and death in the Terminal Complement pathway, doi : https://doi.org/10.1101/2022.04.21.489004 where they have performed CLANS analysis. I understand that CLANS performs blast and visually represents the results as clusters.

What does these clusters mean and is there meaning to the distance between them?

And in the perspective of the paper, what does the CLANS analysis tell the researchers?

Thank you, I am still a student, sorry if this is a simple question. And I have attached a screenshot from the supplementary materials of the article

# Answer
CLANS (Clustering Locus-specific ANnotationS) is a bioinformatics tool that uses BLAST to compare a set of nucleotide or protein sequences to a reference database and visualizes the results as clusters, where each cluster represents a group of similar sequences, the distance between the clusters represents the degree of similarity between the sequences in each cluster. In the context of the paper you are reading, the CLANS analysis likely helped the researchers to identify functional domains within the Terminal Complement pathway, and the clusters formed by the CLANS analysis could represent different functional domains or sub-systems within the pathway, and the distance between the clusters could indicate the degree of similarity or dissimilarity between these domains.
> 1 votes

---
Tags: phylogeny, blast, clustering, sequence-homology
---


# Question
Title: how to download all metadata for a single TCGA dataset and link to data file UUIDs?
I've been able to download specific data types from the TCGA portal (https://portal.gdc.cancer.gov/) using the Manifest files and the GDC Tool.

But it's not clear to me how to download the metadata (e.g. cancer type, stage, etc.) with the links to the UUIDs of a specific data type, e.g. the UUID of a methylation data file to the case UUID:

https://portal.gdc.cancer.gov/files/14dc6f1f-21d0-4b05-9ad5-1d52f5020b12

Any ideas?

# Answer
> 2 votes
To download metadata associated with a specific data file, such as the cancer type and stage, you can use the GDC API (Application Programming Interface). The API allows you to query the GDC database and retrieve information in a structured format, such as JSON or XML.

To download the metadata for a specific data file, you will need the file's unique identifier, also known as its UUID. In your example, the UUID of the methylation data file is "14dc6f1f-21d0-4b05-9ad5-1d52f5020b12".

To use the GDC API, you can construct a query using the GDC API endpoint and the file's UUID. For example, to retrieve the metadata for the methylation data file with the UUID "14dc6f1f-21d0-4b05-9ad5-1d52f5020b12", you could use the following query:

https://api.gdc.cancer.gov/files/14dc6f1f-21d0-4b05-9ad5-1d52f5020b12?expand=metadata

This query uses the "files" endpoint of the GDC API and specifies the file's UUID. The "expand=metadata" parameter tells the API to return the metadata associated with the file. The response will be in json format:

```
{"data": {"data_format": "TXT", "access": "open", "file_name": "5b6353ee-af39-4bb4-a801-1454dfbc338b.methylation_array.sesame.level3betas.txt", "submitter_id": "d67719e6-4149-4492-a265-1f5d5689b904", "data_category": "DNA Methylation", "acl": ["open"], "type": "methylation_beta_value", "platform": "Illumina Human Methylation 450", "file_size": 13168876, "created_datetime": "2022-02-15T17:17:23.357435-06:00", "md5sum": "ab0c90db97a3e12634ef3f54f3ec835d", "updated_datetime": "2022-02-16T19:13:44.087132-06:00", "file_id": "14dc6f1f-21d0-4b05-9ad5-1d52f5020b12", "data_type": "Methylation Beta Value", "state": "released", "experimental_strategy": "Methylation Array", "version": "1", "data_release": "32.0 - 36.0"}, "warnings": {"expand": "unrecognized values: [metadata]"}}

```

---
Tags: public-databases, cancer
---