# Question
Title: From RNA-seq raw counts to normalized data
I have several RNA-seq datasets. Some of them provide RNA-seq raw counts, some provide FPKM, RPKM and some have transcripts per million (TPM) data. Non of them provide fastq files, all data is processed already. At the end I want all datasets to be normalized to TPM.

I'm using this code in order to normalize raw counts to TPM: (using R)

```
rpkm <- apply(X = subset(dataset),
                MARGIN = 2,
                FUN = function(x) {
                  10^9 * x / genelength / sum(as.numeric(x))
                })

TPM <- apply(rpkm, 2, function(x) x / sum(as.numeric(x)) * 10^6) %>% as.data.frame()

```

When RPKM is provided, and no raw counts is available, I use the second line in the same code:

```
TPM <- apply(rpkm, 2, function(x) x / sum(as.numeric(x)) * 10^6) %>% as.data.frame()

```

And when FPKM is provided, I use this formula to transform the data to TPM:

```
TPM = FPKM*X

```

where `X = 1e6/[sum of all FPKM of a sample]`

NOTE: `genelength` is obtained using the `biomart` package in R, to get the transcript length directly out of ensemble.

However, because I know the steps of TPM normalization in theory, one should firstly normalize to gene length, and then to gene depth. I'm not sure this is what the code does, feeling very skiptical about it even thought this code was given to me by a Phd student.

Can you guys please help me with this question? is my code right or should I alter something?

Thank you!

# Answer
> However, because I know the steps of TPM normalization in theory, one should firstly normalize to gene length, and then to gene depth.

No, you should normalize to **transcript length**. Which requires you to know the relative abundance of all the transcripts so you can correct for all their lengths.

I do not think it's possible to un-normalize RPKM or TPM without knowing all the lengths involved in doing the calculations.

Let me give an example: here is RSEM output for just one gene. Note that the software thinks there are three transcripts present for this one gene:

transcript\_id gene\_id length effective\_length expected\_count TPM FPKM IsoPct

```
ENST00000373020 ENSG00000000003 3768    3718.54 179.95  13.88   7.94    48.61

ENST00000494424 ENSG00000000003 820 770.54  0.00    0.00    0.00    0.00

ENST00000496771 ENSG00000000003 1025    975.54  11.94   3.51    2.01    12.29

ENST00000612152 ENSG00000000003 3796    3746.54 0.00    0.00    0.00    0.00

ENST00000614008 ENSG00000000003 900 850.54  33.11   11.17   6.39    39.10

ENST00000373031 ENSG00000000005 1205    1155.54 0.00    0.00    0.00    0.00

ENST00000485971 ENSG00000000005 542 492.54  0.00    0.00    0.00    0.00

```

So, what would you say the effective length should be for the gene level calculations, and would you get that right without having this transcript-level breakdown?
> 2 votes

---
Tags: genetics
---


# Question
Title: Reduce Overfit by removing insignificant proteins through PSEA analysis
I'm training a ML model for disease prediction using protein composition as input, but overfit is present. While looking to *remove proteins and reduce multi-correlation*, as remove composition predictors through PCA have worked but not enough, this problem appeared and managed to use Metascape (1) for it.

I'm using the **Protein Set Enrichment Analysis** (**PSEA**, GSEA for genes) for significance of function as is a method to identify classes of genes or proteins that are over-represented in a large set of genes or proteins, and may have an association with disease phenotypes (2). So:

* *A higher significance coefficient (log10pvalue) between two proteins suggest they have a similar function so one should be removed, or a lower coef shows the need to erase this quasi-independent protein from the analysis?*
* Is this analysis good enough to reduce overfitting by remove redundant observations?

**References:**

(1) https://metascape.org/gp/index.html#/main/step1

(2) https://www.sciencedirect.com/science/article/pii/S2376999819300406

# Answer
> 1 votes
Overfitting is usually resolved by adding more data and not removing data sets. The technical blurb is 'augmenting' your training data set to revolve overfitting.

It is possible there is a very large multi-gene family involved distorting the classification, thus reducing that could reduce the overfit and cause the algorithm to spread the classification over a much broader group of genes. The general advice is to find an alternative genomic source of proteins and add that, not for this specific question but in reducing overfit *per se*.

There are technical issues such as:

* a parameter split, does the training have this? If not its strongly recommended - depending on the algorithm.
* algorithm, for example decision trees are known for overfitting\*, whilst naive Bayes would find difficulty overfitting\**.

I personally think multiple measures would resolve this. However, general advice is to look at alternative methods and training approaches before looking at slimming down a data set. I would suggest there needs to be good evidence - which would need to be formally presented. You could be 100% correct (probably are), but its not an orthodox approach.

Clearly there is information here that isn't presented. What you've done is a feature selection based on genes and we don't know what that result was nor the gene group that emerged from that analysis - so it is difficult to dissect this precisely. It could be there is a compelling biological reason for concern of overfitting based on the result, but we'll not know that. What I am certain about is the orthodoxy in resolving overfitting is not leveraged here and any data scientist would be aware of this.

\*, this is why random forests are considered strong

\**, there are a lot of different approaches to ML and how they are classifying is very different, some are more prone to overfitting than others. Naive Bayesian and 'linear regression' are the least sensitive to overfitting.

---

Please be aware I do ML within the Python framework and its one of the few areas I'm not aware of the ease of cross-over into other platforms. This is simply because its an analysis where Python excels.

---
Tags: machine-learning, proteomics, features
---


# Question
Title: How to modify execution orders for Snakemake?
My goal is to process two samples, `S1` and `S2`, using two rules, `step1`, `step2`. The Snakemake file is like:

```
SAMPLES = ['S1', 'S2']

rule all:
  input: expand("{sample}.done", sample = SAMPLES)

rule step1:
  input: "{sample}"
  output: "{sample}.step1.done"
  shell: "touch {output}"

rule step2:
  input: "{sample}.step1.done"
  output: "{sample}.done"
  shell: "touch {output}"

```

The current execution order is:

```
step1 on S1
step1 on S2
step2 on S1
step2 on S2

```

Is there a way to change the order to:

```
step1 on S1
step2 on S1
step1 on S2
step2 on S2

```

Basically, I want to complete `S1`, then `S2`. Is there a way to specify the order of executions in Snakemake?

# Answer
> 5 votes
I think you should use the priority directive to give precedence to downstream rule(s).

Here below I give higher priority to `step2` so when the first sample has completed `step1` snakemake will run step2 on that sample rather than submitting another sample to step1.

```
SAMPLES = ['S1', 'S2', 'S3']

wildcard_constraints:
    sample= '|'.join([re.escape(x) for x in SAMPLES])

rule all:
    input: expand("{sample}.done", sample = SAMPLES)

rule step1:
    priority: 1
    output: "{sample}.step1.done"
    shell:
        r"""
        sleep 5
        touch {output}
        """

rule step2:
    priority: 10
    input: "{sample}.step1.done"
    output: "{sample}.done"
    shell:
        r"""
        sleep 5
        touch {output}
        """

```

(Note that without `wildcard_constraints` your/this code throws `AmbiguousRuleException`)

# Answer
> 1 votes
> it seems cumbersome to create them to ensure execution order in my case

It took me a while to figure out, but it did not seem cumbersome to me.

You can re-use your rules, see below, and make Step 1 of S2 depend on Step 2 of S1 like so:

```
SAMPLES = ['S1', 'S2']

wildcard_constraints:
    sample = "S[12]"

rule all:
  input: expand("{sample}.done", sample = SAMPLES)

rule step1:
  input: data="{sample}"
  output: data="{sample}.step1.done"
  shell: "touch {output.data}"

rule step2:
  input: data="{sample}.step1.done"
  output: data="{sample}.done"
  shell: "touch {output.data}"

# step2 of S1 has the flag file as output

use rule step2 as step2ForS1 with:
    input: data="S1.step1.done"
    output: data="S1.done", flag=touch("flags/step2-s1-done")

# step1 of S2 needs the flag file as input

use rule step1 as step1ForS2 with:
    input: data="S2", flag="flags/step2-s1-done"
    output: data="S2.step1.done"

```

# Answer
> 0 votes
The problem that you are encountering, is that when snakemake determines the order to run a set of rules it first builds a Directed Acyclic Graph (DAG) of all the rules. This shows which rules are are dependant on each other, and what order they can be safely executed in, without breaking downstream dependencies.

In general, if you want tasks to execute one after another which do not share an explicit file dependency, you can use a Flag File as an ouput for `rule 1` then reference it as an input in `rule 2`.

Unfortunately, I do not think this will help you in your case. You are expanding on the `SAMPLES` list and creating a wildcard. in the `DAG` this creates two competely separate execution branches which have no dependencies with each other. `step1` and `step2` are guaranteed to occur in order for each sample, however there is no guarantees between samples.

Is there a specific reason why you need all of `S2` tasks to be run after `S1`? If so, then your best bet would be to create explicit rules for sample 1 and sample 2, using the `Flag File` to incidate when `step2` for `S1` is done.

---
Tags: snakemake
---


# Question
Title: Is it possible to integrate multiple gene expression datasets and use it for WGCNA?
I have 8 RNA-seq datasets and am interested in looking at genes co-expressed with a specific gene.

Among 8 RNA-seq datasets, 6 have less than 20 samples. Rather than working on each dataset individually for WGCNA, do you think I can merge all datasets and use them for WGCNA?

# Answer
> 2 votes
Yes, if the expression datasets are from the same platform, you can merge all the datasets for common genes among them. However, if they are from different platforms then some batch correction needs to be done before analysis. You can plot a PCA to check how much variability there is among datasets from different platforms and then decide whether batch correction is needed.

# Answer
> 1 votes
I would try to process all data sets in a similar way, then combine them and batch-correct, e.g., using ComBat or WGCNA's own empiricalBayesLM.

---
Tags: r, rna-seq, wgcna
---


# Question
Title: Correct way of merging scRNAseq datasets of healthy and tumor cell?
While merging datasets from similar biological conditions, but different experimental conditions is a well-studied topic, and there are many batch correction techniques available. What I want to understand is the right way to merge scRNA seq counts data of healthy and tumor cells.

Does individually log-normalising the counts, followed by a simple concatenation of the objects works, or there is a better way to do this?

# Answer
I'm aware that you've tagged this as python, but unfortunately I'm most familiar with doing single cell analysis using R. Perhaps there's a python way to do something similar.

Seurat has recently included a workflow for reciprocal PCA integration, which the authors claim should work well when the combined datasets are more different than usual:

> RPCA-based integration runs significantly faster, and also represents a more conservative approach where cells in different biological states are less likely to ‘align’ after integration. We therefore, recommend RPCA during integrative analysis where:
>
> * A substantial fraction of cells in one dataset have no matching type in the other
> * Datasets originate from the same platform (i.e. multiple lanes of 10x genomics)
> * There are a large number of datasets or cells to integrate (see INSERT LINK for more tips on integrating large datasets)

A tutorial workflow can be found here:

https://satijalab.org/seurat/articles/integration\_rpca.html
> 1 votes

---
Tags: python, scrnaseq
---


# Question
Title: Is it common to get different number of SNVs+Indels across samples from vcf files generated using GATK and DRAGEN (counts are higher for GATK)?
We have 2 vcf (Whole Exome Sequencing (**WES**) data; **germline** samples) files (e.g., **vcf\_1** and **vcf\_2**). **vcf\_1** was generated (Ref. genome: hg38) using the **GATK** pipeline for 250 children and their parents and **vcf\_2** was generated (Ref. genome: hg38) using the **DRAGEN** pipeline for another 50 children and their parents (Samples are different). These are raw vcf files and I have not applied any quality-related false positive filters to them.

Now, when we compare the sample-wise SNV counts from **vcf\_1** and **vcf\_2**, there is a huge difference. **vcf\_2** contains on average 4-6 times lesser number of SNVs per sample than **vcf\_1**.

**N.B:** I found this article (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4418901/), where it has been reported that the mean numbers of single-nucleotide variants (SNVs) and small insertions/deletions (indels) detected per sample was around 98k, respectively, for WES. For our case, in the case of GATK-generated vcf samples, this varies from 130k to 188k, and for DRAGEN-generated vcf samples, this varies from 30k to 40k.

I understand both pipelines (GATK and DRAGEN) implement different methods to generate vcf files but is it common/normal to examine that much difference in SNV counts? I would appreciate it if someone can give some explanation (or provide some materials) about it.

# Answer
> 0 votes
I suspect what's happening here is an unfortunate side-effect of making variant-finding algorithms memory-efficient, so that they can work on extremely large populations. By doing the processor-intensive variant calling separately on each individual, calculations can be easily parallelised to run quickly.

For the numbers you have (i.e. 300 children and their parents), there's an alternative approach for generating variant calling that uses all mapped files as input, rather than only one at a time, i.e. instead of running:

```
bcftools mpileup [options] in1.bam
bcftools mpileup [options] in2.bam
bcftools mpileup [options] ...

```

The all-at-once approach has a single mpileup call:

```
bcftools mpileup [options] in1.bam in2.bam ...

```

This will call multiple samples at the same time, so that the output will have identical variant lines for all individuals. It also has the added benefit of getting more information on variation, allowing more rare variants to be identified.

I'm not sure which approach is being used for DRAGEN vs GATK, but I do know that at some stage in the past the GATK standard practice was calling one sample at a time. Curiously, if this were still the case, I'd expect the results to be the other way round, as the one-at-a-time approach tends to be less sensitive.

I guess it's also possible that the pipelines have different sensitivity / specificity tweaks, that could also explain the differences in variant counts.

---
Tags: vcf, ngs, gatk, snv, wes
---


# Question
Title: How to identify sequence origin in DNA shuffle reads?
**The data**

I have a large number of reads from sequences that were generated by randomly shuffling regions of two parent sequences together. See the following image:

The regions that are shuffled are around 100-200 bp in length. The full shuffled sequence is around 1000 bp long. In general, the regions will be in the same order as the original sequence. The shuffled sequences are fragmented and sequenced into 150bp reads. One read can contain more than two borders between shuffle regions from a different parent sequence.

The sequences are very similar (95%). Because of this high similarity you cannot pinpoint the border between two regions. The border will be a range between a variant that is unique to source1 and a variant that is unique to source2. That is ok because I am not interested in knowing from which source an identical region came.

**The question**

I would like to identify which part of the read originates from parent sequence one or two. What would be a suitable tool to do this?

**My attempts**

I tried to look at some multiple sequence alignment tools like MAFFT or Clustal but I am not interested in aligning all reads to each other. They also don't identify the shuffle regions well in my testing.

# Answer
> 2 votes
This is the sort of thing that LAST has been specifically designed to work well for:

> LAST is designed for moderately large data (e.g. genomes, DNA reads, proteomes). It's especially good at finding rearrangements and recombinations: we believe last-split does that more rigorously than anything else.

Here's my usual approach for mapping with LAST:

```
lastdb -uRY4 -R01 reference_sequences.fa reference_sequences.fa
# training is probably not needed if your sub-sequences have no
# variation from the parental sequences
last-train reference_sequences.fa input_sequences.fa > trained.mat
lastal -P 10 -p trained.mat --split reference_sequences.fa input_sequences.fa > aligned.maf

```

LAST is *very* customisable. If you want to adjust the mapping or split sensitivity, have a look at the `lastal` documentation here, and the `last-split` documentation here.

For example, there's a mismap probability adjustment which excludes alignments with a mismap probability above a certain threshold:

```
# only form alignments if the mismap probability is 0.1% or less
lastal --split-m 0.001 ...

```

---
Tags: sequence-alignment, genomics, multiple-sequence-alignment
---


# Question
Title: Get genome coordinates from UniProt entry?
I'm trying to get genome coordinate information programmatically for multiple genes.

As an example, I'd like to get the annotation (list of exons and their chromosomes, starts, and ends) for this UniProt entry: https://www.uniprot.org/uniprotkb/B6TU39/entry.

There is a corresponding GenBank entry, which gives the CDS (https://www.ncbi.nlm.nih.gov/nuccore/EU968504) but I can't see anywhere how to extract the genome coordinates. How can I do this?

Thanks!

# Answer
> 2 votes
You don't want the EU968504, that's the mRNA not the genome sequence. You want https://www.ncbi.nlm.nih.gov/nuccore/CM007647, also linked to from the UniProt page, just under the link to EU968504:

You also need the gene name, `ZEAMMB73_Zm00001d029604` in this case, available further up on the page in the "Names & Taxonomy" section:

Next, we can download the GenBank file for the relevant chromosomal contig using NCBI's entrez API:

```
curl "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=CM007647&rettype=gbwithparts&retmode=text" > CM007647

```

This will take a few minutes. Once that is done, you can parse the file and extract the start and end coordinates using something like this:

```
#!/usr/bin/env python3
from Bio import SeqIO
import sys

gb_record=SeqIO.read(open(sys.argv[1],"r"), "genbank")

for feature in gb_record.features:
    if feature.type == "gene" and sys.argv[2] in feature.qualifiers['locus_tag']:
        if feature.location.strand == 1:
            strand = "+"
        elif feature.location.strand == -1:
            strand = "-"
        else:
            print("Unknown strand: ", feature.location.strand)
            sys.exit(1)

        print("%s\t%s\t%s\t%s\t%s\t%s\t%s\t.\t." % (gb_record.id, "GenBank",
                                                    feature.type, feature.location.start,
                                                    feature.location.end, ".", strand))

```

You pass the GenBank file and gene name as arguments, and this is the result:

```
$ foo.py CM007647  ZEAMMB73_Zm00001d029604
CM007647.1  GenBank gene    78770080    78773320    .   +   .   .

```

# Answer
> 0 votes
You could align back the peptide sequence to the genome, respecting splicing, for example using miniprot from Heng Li.

---
Tags: uniprot, genbank, sequence-annotation
---


# Question
Title: Retrieval protein definition from accession ID
I have a list of protein accessions:

```
XP_025107904.1
PVD22034.1
AJA37836.1

```

How can I obtain the definition of the protein, and the related organism?

Output:

```
XP_025107904.1  receptor-type tyrosine-protein phosphatase kappa-like isoform X2  Pomacea canaliculata

```

Thanks for the help,

Marco

# Answer
> 1 votes
You can do this using the Entrez API. For example, using the three identifiers from your example, you can do:

```
curl "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=protein&id=XP_025107904.1,PVD22034.1,AJA37836.1" > file.xml

```

This will bring an XML file with the following contents:

```
$ cat file.xml
<?xml version="1.0" encoding="UTF-8" ?>
<!DOCTYPE eSummaryResult PUBLIC "-//NLM//DTD esummary v1 20041029//EN" "https://eutils.ncbi.nlm.nih.gov/eutils/dtd/20041029/esummary-v1.dtd">
<eSummaryResult>
<DocSum>
    <Id>1397721498</Id>
    <Item Name="Caption" Type="String">XP_025107904</Item>
    <Item Name="Title" Type="String">receptor-type tyrosine-protein phosphatase kappa-like isoform X2 [Pomacea canaliculata]</Item>
    <Item Name="Extra" Type="String">gi|1397721498|ref|XP_025107904.1|[1397721498]</Item>
    <Item Name="Gi" Type="Integer">1397721498</Item>
    <Item Name="CreateDate" Type="String">2018/06/07</Item>
    <Item Name="UpdateDate" Type="String">2018/06/07</Item>
    <Item Name="Flags" Type="Integer">512</Item>
    <Item Name="TaxId" Type="Integer">400727</Item>
    <Item Name="Length" Type="Integer">1271</Item>
    <Item Name="Status" Type="String">live</Item>
    <Item Name="ReplacedBy" Type="String"></Item>
    <Item Name="Comment" Type="String"><![CDATA[  ]]></Item>
    <Item Name="AccessionVersion" Type="String">XP_025107904.1</Item>
</DocSum>

<DocSum>
    <Id>1383748625</Id>
    <Item Name="Caption" Type="String">PVD22034</Item>
    <Item Name="Title" Type="String">hypothetical protein C0Q70_17837 [Pomacea canaliculata]</Item>
    <Item Name="Extra" Type="String">gi|1383748625|gb|PVD22034.1||gnl|WGS:PZQS|C0Q70_17837[1383748625]</Item>
    <Item Name="Gi" Type="Integer">1383748625</Item>
    <Item Name="CreateDate" Type="String">2018/04/26</Item>
    <Item Name="UpdateDate" Type="String">2018/04/26</Item>
    <Item Name="Flags" Type="Integer">0</Item>
    <Item Name="TaxId" Type="Integer">400727</Item>
    <Item Name="Length" Type="Integer">602</Item>
    <Item Name="Status" Type="String">live</Item>
    <Item Name="ReplacedBy" Type="String"></Item>
    <Item Name="Comment" Type="String"><![CDATA[  ]]></Item>
    <Item Name="AccessionVersion" Type="String">PVD22034.1</Item>
</DocSum>

<DocSum>
    <Id>735659541</Id>
    <Item Name="Caption" Type="String">AJA37836</Item>
    <Item Name="Title" Type="String">fucolectin-related molecule [Littorina littorea]</Item>
    <Item Name="Extra" Type="String">gi|735659541|gb|AJA37836.1|[735659541]</Item>
    <Item Name="Gi" Type="Integer">735659541</Item>
    <Item Name="CreateDate" Type="String">2014/12/25</Item>
    <Item Name="UpdateDate" Type="String">2014/12/25</Item>
    <Item Name="Flags" Type="Integer">0</Item>
    <Item Name="TaxId" Type="Integer">31216</Item>
    <Item Name="Length" Type="Integer">722</Item>
    <Item Name="Status" Type="String">live</Item>
    <Item Name="ReplacedBy" Type="String"></Item>
    <Item Name="Comment" Type="String"><![CDATA[  ]]></Item>
    <Item Name="AccessionVersion" Type="String">AJA37836.1</Item>
</DocSum>

</eSummaryResult>

```

You can now parse this in your favorite language/parser and get the information you need. I... am not too good with XML so my attempts with `xmllint` failed, but here's a naive awk approach:

```
$ awk -F'[][><]' '{ if(/Name="Title"/){desc=$3; spc=$4}else if(/Name="AccessionVersion"/){name=$3} if($0=="</DocSum>"){print name,desc,spc; name=""; desc=""; spc=""}}' file.xml
XP_025107904.1 receptor-type tyrosine-protein phosphatase kappa-like isoform X2  Pomacea canaliculata
PVD22034.1 hypothetical protein C0Q70_17837  Pomacea canaliculata
AJA37836.1 fucolectin-related molecule  Littorina littorea

```

I stress that this is a very bad idea and you should use a proper XML parser, but it might work...

---
Tags: proteins, phylogenetics, uniprot, identifiers
---


title: Snakemake MissingRuleException

text:
# Question
Title: Snakemake MissingRuleException
I'm trying to run a snakefile for the first time with limited coding experience using salmon to index a reference genome. I'm not too sure what I'm doing wrong based on this error message.

```
Error: Building DAG of jobs...
MissingRuleException:
No rule to produce snakefile (if you use input functions make sure that they don't raise unexpected exceptions).

```

Code:

```
rule salmon_index:
        input:
                R1="sequences/*R1.fastq.gz",
                R2="sequences/*R2.fastq.gz"
        output:
                Index="results/salmon_index"
        shell:
                """
                salmon index --threads 8 -t references/genome.fa.gz -i salmon_index
                """

```

# Answer
You're not using wildcards the way snakemake is equipped to work with them. See:

1. FAQs such as https://snakemake.readthedocs.io/en/stable/project\_info/faq.html#how-do-i-run-my-rule-on-all-files-of-a-certain-directory and https://snakemake.readthedocs.io/en/stable/project\_info/faq.html#i-don-t-want-expand-to-use-the-product-of-every-wildcard-what-can-i-do
2. https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#wildcards
> 2 votes

---
Tags: rna-seq, python, snakemake, conda, salmon
---


# Question
Title: subset from a dataframe using other dataframe in R
I have a data frame (named as df) with single column of protein IDs and other file with two columns of protein IDs (named as df1). I need to find the df IDs present in either first or second column of df1 and extract corresponding Id from the df1.

```
df = data.frame(p = c("ENSG00000198400","ENSG00000066044","ENSG00000188612","ENSG00000177606","ENSG00000146648"))
df1 = data.frame(p1 = c("ENSG00000105647","ENSG00000110395","ENSG00000165119","ENSG00000160293","ENSG00000198400","ENSG00000198400"),p2 = c("ENSG00000198400","ENSG00000198400","ENSG00000198400","ENSG00000142949","ENSG00000112936","ENSG00000162434"))

```

Sample output:

```
df2= data.frame(p = c("ENSG00000198400","ENSG00000198400","ENSG00000198400","ENSG00000198400","ENSG00000198400"),p2 = c("ENSG00000105647","ENSG00000110395","ENSG00000165119","ENSG00000112936","ENSG00000162434"))

```

# Answer
you can try this:

```
df = data.frame(p = c("ENSG00000198400","ENSG00000066044","ENSG00000188612",
                       "ENSG00000177606","ENSG00000146648"))
df1 = data.frame(p1 = c("ENSG00000105647","ENSG00000110395","ENSG00000165119",
                         "ENSG00000160293","ENSG00000198400","ENSG00000198400"),
                 p2 = c("ENSG00000198400","ENSG00000198400","ENSG00000198400",
                         "ENSG00000142949","ENSG00000112936","ENSG00000162434"))

df2 = rbind(df1[df1$p1 %in% df$p,],df1[df1$p2 %in% df$p,])

df2 <- transform(df2, p2 = ifelse(df2$p2%in%df$p, p1, p2), p1 = ifelse(df2$p2%in%df$p, p2, p1))

```
> 2 votes

---
Tags: r, dataframe
---


# Question
Title: Error in dimnamesGets(x, value) when trying to read data using Seurat package
*This question has also been asked on Biostars*

I am trying to create a Seurat object using the package "Seurat". When I am reading my raw files using the function `Read10X`, I am getting this error

```
Error in dimnamesGets(x, value) : length of Dimnames[[1]] (45948) is not equal to Dim[1] (45947).

```

I am unsure how to resolve this issue.

R sessionInfo:

R version 4.1.2 (2021-11-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] sp_1.5-0           SeuratObject_4.1.1 Seurat_4.1.1       Matrix_1.5-1

loaded via a namespace (and not attached):
  [1] nlme_3.1-153          spatstat.sparse_2.1-1 matrixStats_0.62.0    RcppAnnoy_0.0.19
  [5] RColorBrewer_1.1-3    httr_1.4.4            sctransform_0.3.4     tools_4.1.2
  [9] utf8_1.2.2            R6_2.5.1              irlba_2.3.5           rpart_4.1-15
 [13] KernSmooth_2.23-20    uwot_0.1.14           mgcv_1.8-38           rgeos_0.5-9
 [17] DBI_1.1.3             lazyeval_0.2.2        colorspace_2.0-3      tidyselect_1.1.2
 [21] gridExtra_2.3         compiler_4.1.2        progressr_0.11.0      cli_3.4.0
 [25] plotly_4.10.0         scales_1.2.1          lmtest_0.9-40         spatstat.data_2.2-0
 [29] ggridges_0.5.3        pbapply_1.5-0         goftest_1.2-3         stringr_1.4.1
 [33] digest_0.6.29         spatstat.utils_2.3-1  pkgconfig_2.0.3       htmltools_0.5.3
 [37] parallelly_1.32.1     fastmap_1.1.0         htmlwidgets_1.5.4     rlang_1.0.5
 [41] rstudioapi_0.14       shiny_1.7.2           generics_0.1.3        zoo_1.8-10
 [45] jsonlite_1.8.0        spatstat.random_2.2-0 ica_1.0-3             dplyr_1.0.10
 [49] magrittr_2.0.3        patchwork_1.1.2       Rcpp_1.0.9            munsell_0.5.0
 [53] fansi_1.0.3           abind_1.4-5           reticulate_1.26       lifecycle_1.0.2
 [57] stringi_1.7.6         MASS_7.3-54           Rtsne_0.16            plyr_1.8.7
 [61] grid_4.1.2            parallel_4.1.2        listenv_0.8.0         promises_1.2.0.1
 [65] ggrepel_0.9.1         deldir_1.0-6          miniUI_0.1.1.1        lattice_0.20-45
 [69] cowplot_1.1.1         splines_4.1.2         tensor_1.5            pillar_1.8.1
 [73] igraph_1.3.4          spatstat.geom_2.4-0   future.apply_1.9.1    reshape2_1.4.4
 [77] codetools_0.2-18      leiden_0.4.3          glue_1.6.2            data.table_1.14.2
 [81] png_0.1-7             vctrs_0.4.1           httpuv_1.6.6          polyclip_1.10-0
 [85] gtable_0.3.1          RANN_2.6.1            purrr_0.3.4           spatstat.core_2.4-4
 [89] tidyr_1.2.1           scattermore_0.8       future_1.28.0         assertthat_0.2.1
 [93] ggplot2_3.3.6         mime_0.12             xtable_1.8-4          later_1.3.0
 [97] survival_3.2-13       viridisLite_0.4.1     tibble_3.1.8          cluster_2.1.2
[101] globals_0.16.1        fitdistrplus_1.1-8    ellipsis_0.3.2        ROCR_1.0-11
```

# Answer
> 1 votes
Check the length of your files. The error clearly says that there is a difference in the length of your files. Your file length should be equal.

I think the problem may be the presence of a header in one of your files. Removing the header may solve the problem.

# Answer
> 1 votes
There is a mismatch of 1 line between the several files, suggesting that probably either the barcodes of features file has a header line which is not standard. Remove any header lines from these files.

---
Tags: r, scrnaseq, seurat, bioconductor, single-cell
---


title: Comparison of fastq files reads

text:
# Question
Title: Comparison of fastq files reads
My goal is to compare reads from two different fastq files on a Linux machine. The following are the comparisons to perform:

1. How many common reads are between the two fastq files?
2. How many reads are present in one fastq file that are not present in the other? i.e reads that are unique to each file.
3. Where and how do the sequences differ in both files? Is there a difference in nucleotide base at the end of the sequence(3' end), the start of the sequence(5'end) or is it somewhere in the middle of the sequence? Is there a difference in length of the sequence? If so, what is the length and what is the difference? I would like to get the original sequences plus the difference between them for analysis.
4. Does the quality score differ? Where and how do the quality of sequences differs?

Any guidance regarding how to approach this question is very appreciated.

Edit: I'm working on paired-end reads from the same individual. What I want to do is to compare different trimming tools from their output.

# Answer
There are several things you can try. For example, Simka (https://github.com/GATB/simka) is tool that compares the profiles of the genomic content of fastq files and generates several types of plots e.g. heatmap.

Another tool that might be useful is multifastqc (https://multiqc.info/), which aggregates results and has multiple options.
> 1 votes

---
Tags: python, fastq, reads
---


# Question
Title: Is there an equivalent of Entrez queries for local BLAST?
I have been having some trouble with running command line BLASTn remotely. Therefore, upon recommendations, and considering I have \>100 searches planned, I am considering just downloading the nt database (thankfully possible).

My issue is that one of the settings that proved valuable with my data in WebBlast is running an Entrez query. I get more useful results for each query when instead of one search, I run three, with each of the sub-searches delimited by subject sequence length (e.g. `-entrez_query 700:1000[slen]` to only get the results which are between 700 and 1000 nucleotides long).

WebBlast and remote Blast have that, local Blast does not.

Since I only really need 2 Entrez queries, I guess it would not be conceptually impossible to get my copy of nt database as FASTA, and then run a script to split it into three directories (700-1000nt, 1000-1400nt, the rest while discarding all\<700 nt), then make three blastable databases out of it. However it does feel rather inelegant. I was wondering if anyone might have a better solution!

# Answer
I looked back at the question and have a more comprehensive solution.

1. There is no option for `-entrez_query` on a local blast search. The information is here

> entrez\_query Restrict search with the given Entrez query. Remote searches only.

There is an option without `-outfmt` for `length` and that is worth investigating.

> length means Alignment length

2. ### Perl one liner solution to NCBI's many parameter problem

a) Preprocess your fasta file into separate databases of the required length ... by turning a fasta flatfile db into one where all the subject sequences are 700 - 1000.

```
perl -0777 -nle 'while (/^(>.+\n^[AGCT-NRYMKWS]{700,1000})\n/gm){print $1,"\n"}' example.fa > example700-1000.fa

```

b) Now make the database on the restricted length fasta flatfile db,

```
makeblastdb -in example700-1000.fa -parse_seqids -blastdb_version 5  -out laura700-1000.db -title "Laura 700 - 1000" -dbtype nucl

```

This db will now contain an equivalent of the `-entrez_query 700:1000[slen]`. Obviously, you would need to make a new db for each change in the sequence length query. Note the option `-taxid_map`.

c) Finally, run the blast search on the restricted length blast db.

```
blastn -query laura.fas -db laura700-1000.db -evalue 1e-6 -num_threads 4 -out blasthits.txt

```

---

#### Code testing

**Input is test.fa**

```
>here
AGAGACGTAGAGAAAAAAAAAA
>here2
TAGAGA
>here3
AGAGACGTAGAGA

```

Modified perl oneliner because sequence length is different,

```
perl -0777 -nle 'while (/^(>.+\n^[AGC-T]{6,13})\n/gm){print $1,"\n"}' test.fa

```

**Output**

```
>here2
TAGAGA
>here3
AGAGACGTAGAGA

```

Thus the sequence that was outside the range 6 to 13 (sequence `here2` and `here3`) is excluded (`here`). This seems alot easier to me that all those `blastn` options.

**Note** NRYMKWS refer to the degenerate nucleotide code
> 1 votes

# Answer
"Often, one needs to search ... a specific subset of sequences within an existing database. For these type of searches a convenient way to conduct them is by creating a virtual BLAST database"

From "Use blastdb\_aliastool to manage the BLAST databases" https://www.ncbi.nlm.nih.gov/books/NBK569848/
> 0 votes

---
Tags: blast, bash, blastn, entrez
---


title: Is there public RESTful api for Gnomad?

text:
# Question
Title: Is there public RESTful api for Gnomad?
I currently find Harvard's RESTful API for ExAC extremely useful and I was hoping that a similar resource is available for Gnomad?

Does anyone know of a public access API for Gnomad or possibly any plans to integrate Gnomad into the Harvard API?

# Answer
> 15 votes
As far as I know, no but the vcf.gz files are behind a http server that supports *Byte-Range*, so you can use `tabix` or any related API:

```
$ tabix "https://storage.googleapis.com/gnomad-public/release-170228/vcf/exomes/gnomad.exomes.r2.0.1.sites.vcf.gz" "22:17265182-17265182"
22  17265182    .   A   T   762.04  PASS    AC=1;AF=4.78057e-06;AN=209180;BaseQRankSum=-4.59400e+00;ClippingRankSum=2.18000e+00;DP=4906893;FS=1.00270e+01;InbreedingCoeff=4.40000e-03;MQ=3.15200e+01;MQRankSum=1.40000e+00;QD=1.31400e+01;ReadPosRankSum=2.23000e-01;SOR=9.90000e-02;VQSLOD=-5.12800e+00;VQSR_culprit=MQ;GQ_HIST_ALT=0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|1;DP_HIST_ALT=0|0|0|0|0|0|0|0|0|0|0|1|0|0|0|0|0|0|0|0;AB_HIST_ALT=0|0|0|0|0|0|0|0|0|0|0|0|1|0|0|0|0|0|0|0;GQ_HIST_ALL=1591|589|120|301|650|589|1854|2745|1815|4297|5061|2921|10164|1008|6489|1560|7017|457|6143|52950;DP_HIST_ALL=2249|1418|6081|11707|16538|9514|28624|23829|7391|853|95|19|1|0|0|1|0|1|0|0;AB_HIST_ALL=0|0|0|0|0|0|0|0|0|0|0|0|1|0|0|0|0|0|0|0;AC_AFR=0;AC_AMR=0;AC_ASJ=0;AC_EAS=0;AC_FIN=1;AC_NFE=0;AC_OTH=0;AC_SAS=0;AC_Male=1;AC_Female=0;AN_AFR=11994;AN_AMR=31324;AN_ASJ=7806;AN_EAS=13112;AN_FIN=20076;AN_NFE=94516;AN_OTH=4656;AN_SAS=25696;AN_Male=114366;AN_Female=94814;AF_AFR=0.00000e+00;AF_AMR=0.00000e+00;AF_ASJ=0.00000e+00;AF_EAS=0.00000e+00;AF_FIN=4.98107e-05;AF_NFE=0.00000e+00;AF_OTH=0.00000e+00;AF_SAS=0.00000e+00;AF_Male=8.74386e-06;AF_Female=0.00000e+00;GC_AFR=5997,0,0;GC_AMR=15662,0,0;GC_ASJ=3903,0,0;GC_EAS=6556,0,0;GC_FIN=10037,1,0;GC_NFE=47258,0,0;GC_OTH=2328,0,0;GC_SAS=12848,0,0;GC_Male=57182,1,0;GC_Female=47407,0,0;AC_raw=1;AN_raw=216642;AF_raw=4.61591e-06;GC_raw=108320,1,0;GC=104589,1,0;Hom_AFR=0;Hom_AMR=0;Hom_ASJ=0;Hom_EAS=0;Hom_FIN=0;Hom_NFE=0;Hom_OTH=0;Hom_SAS=0;Hom_Male=0;Hom_Female=0;Hom_raw=0;Hom=0;POPMAX=FIN;AC_POPMAX=1;AN_POPMAX=20076;AF_POPMAX=4.98107e-05;DP_MEDIAN=58;DREF_MEDIAN=5.01187e-84;GQ_MEDIAN=99;AB_MEDIAN=6.03448e-01;AS_RF=9.18451e-01;AS_FilterStatus=PASS;CSQ=T|missense_variant|MODERATE|XKR3|ENSG00000172967|Transcript|ENST00000331428|protein_coding|4/4||ENST00000331428.5:c.707T>A|ENSP00000331704.5:p.Phe236Tyr|810|707|236|F/Y|tTc/tAc||1||-1||SNV|1|HGNC|28778|YES|||CCDS42975.1|ENSP00000331704|Q5GH77||UPI000013EFAE||deleterious(0)|benign(0.055)|hmmpanther:PTHR14297&hmmpanther:PTHR14297:SF7&Pfam_domain:PF09815||||||||||||||||||||||||||||||,T|regulatory_region_variant|MODIFIER|||RegulatoryFeature|ENSR00000672806|TF_binding_site|||||||||||1||||SNV|1||||||||||||||||||||||||||||||||||||||||||||,T|regulatory_region_variant|MODIFIER|||RegulatoryFeature|ENSR00001729562|CTCF_binding_site|||||||||||1||||SNV|1||||||||||||||||||||||||||||||||||||||||||||

```

**UPDATE**: 2019: the current server for gnomad doesn't support Byte-Range requests.

# Answer
> 6 votes
The new gnomAD site (as of August 2019) says no, no API yet:

> How do I query a batch of variants? Do you have an API?
>
> We currently do not have a way to submit batch queries on the browser, but we are actively working on developing an API for ExAC/gnomAD. If you would like to learn about GraphQL, which we will use to work with the API, an overview can be found at https://graphql.org. You can also obtain information on all variants from the VCFs and Hail Tables available on our downloads page.

But, the web interface itself already makes POST requests to https://gnomad.broadinstitute.org/api to send and receive JSON/GraphQL. So, you can make those same queries programmatically right now, even if it's not officially a public API.

Here's an example in Python to get some basic info on variants for a particular gene. This way you get simple nested Python objects to work with:

```
  { 'consequence': 'intron_variant',
    'pos': 6928442,
    'rsid': 'rs782435448',
    'variant_id': '12-6928442-C-A'},
  { 'consequence': 'splice_region_variant',
    'pos': 6928462,
    'rsid': None,
    'variant_id': '12-6928462-C-A'},
  { 'consequence': 'splice_acceptor_variant',
    'pos': 6928464,
    'rsid': 'rs782577109',
    'variant_id': '12-6928464-G-A'},
  { 'consequence': 'missense_variant',
    'pos': 6928466,
    'rsid': 'rs782208003',
    'variant_id': '12-6928466-C-T'},

```

(I found it useful to go this route because then the full metadata visible in the gnomAD web interface is then available, including the per-variant details like allele counts by population. I couldn't find this information in the other APIs described here.)

# Answer
> 4 votes
You can browse gnomAD variants with ClinGen Allele Registry (there is API spec available).

# Answer
> 2 votes
I faced same issue recently, I found those link and python script:

gnomAD GraphQL api https://gnomad.broadinstitute.org/api It works great but it is a kind of different query language. Please check here for the docs: https://graphql.org/learn/queries/

gnomAD Python Api https://github.com/furkanmtorun/gnomad\_python\_api

# Answer
> 2 votes
I found Jesse's code quite usefull ! For those who try to reproduce it, you should now add the reference genome ID, such as :

```
#!/usr/bin/env python

import requests
# import pprint
# prettyprint = pprint.PrettyPrinter(indent=2).pprint

def fetch(jsondata, url="https://gnomad.broadinstitute.org/api"):
    # The server gives a generic error message if the content type isn't
    # explicitly set
    headers = {"Content-Type": "application/json"}
    response = requests.post(url, json=jsondata, headers=headers)
    json = response.json()
    if "errors" in json:
        raise Exception(str(json["errors"]))
    return json

def get_variant_list(gene_id, dataset="gnomad_r2_1"):
    # Note that this is GraphQL, not JSON.
    fmt_graphql = """
    {
        gene(gene_id: "%s", reference_genome: GRCh38) {
          variants(dataset: %s) {
            consequence
            pos
            rsid
            variant_id: variantId
          }
        }
      }
    """
    # This part will be JSON encoded, but with the GraphQL part left as a
    # glob of text.
    req_variantlist = {
        "query": fmt_graphql % (gene_id, dataset),
        "variables": {"withFriends": False}
        }
    response = fetch(req_variantlist)
    return response["data"]["gene"]["variants"]

prettyprint(get_variant_list("ENSG00000010610"))

```

# Answer
> 1 votes
I created a python package based on SQLite databases, where you can easily query all gnomAD variants for GRCh37/38. https://github.com/KalinNonchev/gnomAD\_DB I have precomputed SQLite databases for gnomAD WGS for GRCh37/38 in the description of the package. Please take a look there."

# Answer
> -1 votes
https://github.com/KalinNonchev/gnomAD\_MAF would solve your problem I think. It is in python and you can annotate fast variants with their allele frequencies.

---
Tags: variation, api
---


title: PL and QUAL values on VCF file?

text:
# Question
Title: PL and QUAL values on VCF file?
I want to filter my VCF file to include only the relevant information I need but I have some questions about the results I got.

This is what the first 10 lines of my VCF file looks like:

```
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  DATA_VARIANT
chrI    826 .   T   C   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    1/1:344,30,0:57:10:0,10,0,0:0,5
chrI    831 .   C   T   71  .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    0/1:139,36,274:71:12:0,8,0,4:3,2
chrI    836 .   G   T   59  .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    0/1:140,49,373:59:16:0,0,12,4:4,1
chrI    848 .   C   T   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    1/1:556,48,0:75:16:0,0,0,16:0,5
chrI    855 .   G   A   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    0/1:399,48,139:59:16:12,0,4,0:1,4
chrI    858 .   ATT AT  255 .   CNV=1;TYPE=INDEL    GT:PL:GQ:DP:ADP:ACN 0/1:480,48,160:85:16:4,12:1,4
chrI    876 .   T   C   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    1/1:865,75,0:102:25:0,25,0,0:0,5
chrI    904 .   G   A   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    1/1:1261,112,0:138:37:37,0,0,0:0,5
chrI    920 .   A   T   67  .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    0/1:244,145,1426:67:48:41,0,0,7:4,1
chrI    925 .   T   C   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    1/1:1786,157,1:255:52:0,52,0,0:0,5

```

**QUAL:** As I understand it, the QUAL is based off of Phred quality scores, which would mean that a QUAL value of anything greater than 20 would be fairly sufficient in terms of quality for that genotype. However, most of my QUAL scores are 255 and from what I've read that means there is no mapping quality available?

* What does this mean for my data? I used the same tool (NGSEP) to map and call the variants on this strain so it should have mapping quality scores.
* Do I need to filter out all the 255 scores? Or can they still be used?

**PL:** From what I've read about PL scores, they indicate the likelihood of the presence of a genotype (0/0, 0/1, 1/1). Based on the results, there should be a clear least likely genotype, "almost" likely genotype, and a most likely genotype indicated by a value of 0.

* Some of my PL values do not include a 0 (i.e., the second variant call = 139,36,274). What does this mean for my data?

Any help would be great! Thank you!

**Note:** I used NGSEP as my variant caller

# Answer
The VCF specification provides the definition for the QUAL field.

However, QUAL values are often capped by variant callers to a given value. I have seen 100 being used and according to this documentation on NGSEP it seems 255 has been chosen for this caller as the maximum value. So you are all good keeping these high QUAL variant sites and filtering only those below a certain value (e.g. 20)

PL is also defined in the VCF specification: phred-scaled genotype likelihoods rounded to the closest integer (and otherwise defined precisely as the GL field).

Again, different variant callers may have different implementations for the specification. GATK for example normalizes the values so that the most likely genotype has a PL of 0. But this might not be followed in other callers. It seems that NGSEP simply does not apply this normalization or does it differently. I think this is a detail you can ask directly the authors here by opening an "issue".

My suggestion is that you rather use GQ (genotype quality) if you want to hard-filter your variants. This is more standard to do, as opposed to filtering using PLs. These values are also in phred-scale.
> 5 votes

---
Tags: vcf, filtering, quality-control, variant-filtration, phred-scores
---


# Question
Title: Allele effects by parental-origin estimated from genotyping array data
We are dealing with relatively large family data (\> 1,000 trios) with genotyping array data. We would like to see the allele effects on children separately of mothers, fathers, and children, as in a previous study (PMID: 34282336). Unfortunately, this is difficult for us to replicate, so we need another method. Is there any software that can achieve this or a similar traditional approach? Also, can that be applied to grandparents?

In my understanding of the above mentioned study is they:

1. identified the parental genomes within the child genome;
2. thereon conducted a GWAS for birth weight

Overall, the aim of the previous study was to estimate the alleles responsible for birth weight or in our situation a comparable phenotype. It is important to note the prior study did not use publicly available algorithms.

Thus, we wish to identify the maternal and paternal alleles within their progeny via public domain algorithms, and therefore do not necessarily need to exactly reproduce the previous study.

We already have imputed data, and are using GCTA for usual GWAS. To achieve our goal, we need to identify the transmitted and non-transmitted alleles in parents.

# Answer
> 1 votes
The first haplotyping step is typically done in two stages:

1. Imputation of missing genotypes - can be done using the Michigan Imputation Server
2. Phasing of called genotypes to create semi-contiguous haplotypes - can be done using SHAPEIT

I have many thoughts about the best approach for GWAS (mostly around the idea of population subsampling), which could potentially also apply to what you want to do. Unfortunately, it's been a long time since I did family population genetics, and I've forgotten the things we did to deal with complex population structure. I'm pretty sure we used PLINK for at least some of it, but I notice that PLINK has a disclaimer about close family structure, which suggests it may not work so well for trio data:

> This method does not properly adjust for small-scale family structure. As a consequence, it is usually necessary to prune close relations with e.g. --king-cutoff before using --glm for genome-wide association analysis. (Note that biobank data usually comes with a relationship-pruned sample ID list; you can use --keep on that list, instead of performing your own expensive --king-cutoff run.) If this throws out more samples than you'd like, consider using mixed model association software such as SAIGE, BOLT-LMM, GCTA, or FaST-LMM instead; or regenie's whole genome regression.

---
Tags: genome, gwas, software-recommendation, microarray
---


# Question
Title: Filtering file with AWK and writing output to new file
Hi I'm trying to filter a column in a PAF file with the AWK command and write it to a new file:

`awk '$12>=20 {print $0}' /PATH/FILE.paf > /OUT/PATH/FITLERED_FILE.paf`

This command makes sure the 12th column has a value of at least 20 and then prints the entire contents of the file - with the filtered column - to a new file.

This works for one file but I have 90+ I need to filter through. I tried using a for loop:

`for i in /PATH/TO/FILES/*_mapped.paf; do awk '$12>=20 {print $0 > "awk_{}.out"}' $i; done`

Which didn't work. And then I tried using a pipe command:

`cat /Users/rimo/Desktop/Pat_Assemblies/flye_mapped.txt | parallel "awk '{$12>=20 print $0}' /Users/rimo/Desktop/Pat_Assemblies/Assemblies_Flye/{} > awk_{}"`

Which also didn't work. I think the problem I am having is format-based but I would appreciate it if someone could steer me in the right direction!

# Answer
> 3 votes
The `{}` isn't in any way special to the shell or `awk`. It is used as a replacement by certain other programs, off the top of my head I can only think of `find` and `xargs` where it is used by convention, but it isn't a shell thing or an awk thing. Here, you are using it inside the `awk` command so it just breaks the awk's format because `{ }` introduces a command block and this one is empty.

Your file name is `$i`, so if you want to make `foo_mapped.paf` into `awk_foo_mapped.paf.out`, you need `awk_$i.out`. Something like this:

```
for file in /PATH/TO/FILES/*_mapped.paf; do
    fileName="${file##*/}"
    awk '$12>=20' "$file" > awk_"$fileName".out
done

```

Note how I removed the `print $0`. This is because the default action for awk when something evaluates to `true` is to print the current line, so `awk '$12>=20 {print $0 }` is equivalent to `awk '$12>=20`.

You can also use your original approach of having `awk` deal with the file name by using awk's internal `FILENAME` variable. Like this:

```
awk '$12>=20{ print > "awk_"FILENAME".out" }' /PATH/TO/FILES/*_mapped.paf

```

Depending on the number of files, you might need to close the file handle after writing. If you get an error about too many open files, try:

```
mawk '$12>=20{print > "awk_"FILENAME".out"; close("awk_"FILENAME".out")}' /PATH/TO/FILES/*_mapped.paf

```

---
Tags: filtering, awk, parallel, command-line, paf
---


# Question
Title: Add file name to FASTA file header using SeqIO
I have several fasta files that represent the proteins from different orthogroups. However, the orthogroup that each of the proteins belongs to is not indicated within the files. I am leveraging BioPython and I want to add this information to the header of each sequence in the fasta files.

```
import os
from pathlib import Path
from Bio import SeqIO

with open('output.fasta', 'w') as output:
    for file in os.listdir(dirname):
        seq = SeqIO.parse(file, 'fasta')
        ortog = Path(file).stem
        for record in seq:
            record.description = ortog
        SeqIO.write(record, output, 'fasta')

```

However, the output always appears with the last sequence duplicated. In the case of three fasta files, the output appears as:

```
>BBC_8_672 OG0015396-PROT
MKFSSVTAITLATVATVATAKKGEHDFTTTLTLSSDGSLTTTTSTHTTHKYGKFNKTSKS
KTPNHTGTHKYGKFNKTSKSKTP
>CAT_6_3128 OG0020551-PROT
MEKIKEKLNSLKLESESWQEKYEELREQLKELEQSNTEKENEIKSLSAKNEQLDSEVEKL
ESQLSDTKQLAEDSNNLRSNNENYTKKNQDLEQQLEDSEAKLKEAMDKLKEADLNSEQMG
RRIVALEEERDEWEKKCEEFQSKYEEAQKELDEIANSLENL
>ZSA_0_3121 OG0025647-PROT
LLLQYQLHPALLPKALLLQYQLHQALLLKALLLQ*PALLLKALLLQYQLHQALLLKALLL
QYQLHPALLLKALLLQYQLHP
>ZSA_0_3121 OG0025647-PROT
LLLQYQLHPALLPKALLLQYQLHQALLLKALLLQ*PALLLKALLLQYQLHQALLLKALLL
QYQLHPALLLKALLLQYQLHP

```

Could anyone tell me why this is happening?

# Answer
> 0 votes
This is because you are also parsing your output file. Remember, the first rule of debugging is to print all the things! Most bugs are because some variable doesn't actually have the value you think it has, and printing them is the best way of finding the problem. I reproduced what (I am guessing) is your setup here:

```
$ ls
OG0015396-PROT.fa  OG0020551-PROT.fa  OG0025647-PROT.fa

```

I then ran your script with an extra `print` statement:

```
#!/usr/bin/env python3
import os
from pathlib import Path
from Bio import SeqIO

with open('output.fasta', 'w') as output:
       for file in os.listdir("."):
              print("FF:",file)
              seq = SeqIO.parse(file, 'fasta')
              ortog = Path(file).stem
              for record in seq:
                     record.description = ortog
              SeqIO.write(record, output, 'fasta')

```

Which printed:

```
$ foo.py
FF: OG0015396-PROT.fa
FF: OG0020551-PROT.fa
FF: output.fasta
FF: OG0025647-PROT.fa

```

As you can see, since the `output.fasta` file is created by the `open()` call, by the time you use `listdir`, it's already there and included. This results in duplicated sequences (not always the last one). I'm sure there are more elegant solutions, but a simple workaround is to skip the output file name:

```
#!/usr/bin/env python3
import os
from pathlib import Path
from Bio import SeqIO

with open('output.fasta', 'w') as output:
       for file in os.listdir("."):
              if file == "output.fasta":
                     continue
              seq = SeqIO.parse(file, 'fasta')
              ortog = Path(file).stem
              for record in seq:
                     record.description = ortog
              SeqIO.write(record, output, 'fasta')

```

Of course, this whole thing is a bit overkill for such a simple task and you could just use `sed` and the shell instead:

```
for file in *fa; do sed "1s|$| ${file//.fa}|" "$file"; done > output.fasta

```

---
Tags: python, phylogenetics, orthologues, multi-fasta
---


title: Understanding PAF format

text:
# Question
Title: Understanding PAF format
Does anyone have any experience with PAF files (minimap2 output)? I don’t fully understand what my output means.

I’ve looked at the minimap2 explanation (https://github.com/lh3/miniasm/blob/master/PAF.md) but I want to use the data here to go back and extract regions in the genome that the gene target mapped to. I just don't know what specific columns to look at to accomplish this.

Here is a little snippet of a PAF file I have:

```
contig_25   522460  514241  515712  +   GENE_TARGET 4988    1962    3649    258   1687  60
contig_26   413461  403718  406719  +   GENE_TARGET 4988    17      3649    1762  3633  60
contig_44   24218   16005   17476   +   GENE_TARGET 4988    1962    3649    258   1687  60
contig_52   244411  236199  237670  +   GENE_TARGET 4988    1962    3649    258   1687  60
contig_57   1043583 1025163 1026895 +   GENE_TARGET 4988    3046    4986    1376  2059  60

```

Is this saying that the gene target mapped to coordinates `514241 through 515712` on `contig_25` but only coordinates `1962 through 3649` of the gene actually mapped to that contig? I guess I’m mostly confused on the descriptions for the coordinates and what they actually mean.

I think for my use case I have to use the gene target as the reference sequence as I am trying to extract homologous regions to sequences of interest. I.e., my project involves 90+ assemblies and I would like to extract the regions in the contigs where they map to the gene target. I used minimap2 to map the assemblies to the gene target, used sort and awk to get the best/high quality matches, converted PAF to BED format, and then used seqtk subseq to extract the sequences from the respective assemblies.

This is what I get when I use the gene target as a reference sequence:

```
GENE    4988    17      3649    +   contig_26   413461  403718  406719  1762    3633    60
...
GENE    4988    3046    4986    +   contig_57   1043583 1025163 1026895 1376    2059    60

```

Would I use columns 8 and 9 as the region range to go back and extract the sequences from?

# Answer
As @user172818 mentioned in a comment to your question, it looks like you might have used the query as the reference; they should be switched round.

A minimap2 mapping command is typically represented as follows:

```
minimap2 reference.fa query.fa > results.paf

```

This follows the "indexed reference first" command line order used in the overwhelming majority of mapping tools. In the minimap2 PAF output (which you have linked), the first column represents the *query* sequence name, and the 6th column represents the *target* sequence name (i.e. the opposite order from the command line order, also following with established convention).

Based on the desired mapping operation that you have described (i.e. "I want to use the data here to go back and extract regions in the genome that the gene target mapped to"), your minimap2 command should probably look like this:

```
minimap2 genome_contigs.fa GENE_TARGET.fa > results.paf

```

This will appropriately map the target region to the regions in the genome, and then columns 8 and 9 will represent "regions in the genome that the gene target mapped to."

Minimap2 (and other mappings programs) match query sequences to a reference, with the usual expectation being that the output will contain a single primary mapping per input query sequence (sometimes with secondary mappings to other locations). In other words, they try to find one match to any reference sequence for each query. Changing the order of mapping to target-as-reference means that the program will attempt to find one match to the target sequence(s) for each input reference sequence, which is usually not what is desired.

In any case, the columns in the PAF file correspond to the closest sequence ID to the *left* of them (which is, confusingly, ordered the opposite way round to the command line order):

```
$ minimap2 [genome_contigs.fa] {GENE_TARGET.fa}
{GENE 4988 17   3649} + [contig_26 413461  403718  406719 ] 1762 3633 60
...
{GENE 4988 3046 4986} + [contig_57 1043583 1025163 1026895] 1376 2059 60

$ minimap2 [GENE_TARGET.fa] {genome_contigs.fa}
{contig_26 413461  403718  406719 } + [GENE 4988 17   3649] 1762 3633 60
...
{contig_57 1043583 1025163 1026895} + [GENE 4988 3046 4986] 1376 2059 60

```
> 5 votes

---
Tags: read-mapping, minimap2, paf
---


# Question
Title: Is picking variants belonging to a specific group of genomes from VCF file possible?
GnomAD has publicly available VCF files that have variant data gathered from over 15k individuals. For a university project, I want to randomly select 1k of these individuals and get the variant data belgoning to them only.

I am very new to bioinformatics in general and upon checking VCF files and how they work, my current conclusion is that this task is not achievable through filtering the vcf file in any way since variants in these files have no direct connection to the genomes they were found in (no array of individual IDs in a variant entry that show which individuals had this variant).

But this task was given to me by a professor whom I believe is competent to realize the task was impossible if it is like I thought, so is there something I am missing?

How can I go about randomly selecting a 1000 individuals from this 15 000 pool, and then generating a new VCF file with the data from these selected 1000?

Are there any fields that show which individuals genome samples showed a spesific variant entry in the VCF file?

# Answer
It is. A VCF ought to be sub-settable to a specific list of samples. The only thing you have to come up with for your specific case, is a list of 1K random sample IDs within the original VCF.

Geenerally speaking, knowing that you can do this and how involves getting familiar with the **VCF specification** as well as the main bioinfo tool used to manipulate and work with this format: **BCFtools**

Let's use **1000 genomes data** as an example.

You can use a bit of UNIX together with BCFtools in order to get a list of 1K random IDs from individuals:

(sort -R instead of shuf might work too)

```
bcftools query -l ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz | shuf | head -1000 > myRandomIDs.txt

```

Then again BCFtools to subset the list of individuals from the original VCF and into a new VCF file:

```
bcftools view --samples-file myRandomIDs.txt ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz -o myNewVCF.vcf

```

If you dig a bit deeper into the BCFtools documentation, you will find many things you can do extra if you wanted, like re-calculating the allele frequencies based on this new cohort of individuals, if you wanted that too. Or subset for positions/variants, instead of individuals.

**GnomAD data**

Using GnomAD data, selecting 1K random samples is impossible since data is aggregated and instead of genotypes per individual (and individuals you can select thus) we only have the frequency of the variants across various sub-cohorts within the GnomAD project. No individual level data.

The most you can do is select one of such sub-cohorts and filter for variants with a frequency \> 0 (or a range) in such cohort. You could get variants present only and only is a specific cohort too, by using either frequency information or other info such as allele counts, etc.

The example below filters for variants with a frequency above 0 in the South Asian cohort and 0 frequency in the Southern European cohort. You can have a look at the header of the VCF in order to get to know all cohorts included and the meaning of all the information contained within the INFO filed in the VCF

```
bcftools filter -i "AF_sas > 0 & AF_nfe_seu = 0" gnomad.exomes.r2.1.1.sites.21.vcf.bgz

```

But I'm afraid that's not exactly what you wanted to do in the first place.
> 0 votes

---
Tags: vcf, gnomad
---


# Question
Title: How do I generate a variant list using Illumina reads from a Salmonella genome?
I am planning on performing phylogenetic analysis of Salmonella specimens using WGS data from PulseNet and GenomeTrakr, for the purpose of public health surveillance & to provide context for foodborne outbreaks.

I have been using the Nextstrain workflow to generate and visualize trees for SARS-CoV-2 specimens, and I am planning on adapting the bacterial workflow outlined here for my bacterial isolates: https://docs.nextstrain.org/en/latest/tutorials/creating-a-bacterial-pathogen-workflow.html

What approach would you use to process my reads into VCF files for phylogenetic analysis? I've been referencing information from this similar post from 5 years ago: How do I generate a variant list (i.e. VCF file) using Illumina reads from a human genome?

How much of this information is still current/relevant and applicable to Salmonella isolates?

# Answer
I would express caution about cross-over because for both SARS-CoV-2 and the human genome you can obtain a good assembly from a reference genome. The SARS-CoV-2 genome is so small and read depth high enough this is not needed however. SARS-CoV-2 there are critical issues in assemblies for certain regions which have not been addressed, but for CDS there's no real difficulties here.

A reference genome assembly is not good practice for what are often environmental bacteria (bacteria *per se*) and *de novo* assembly is strongly recommended, even if it is to complement a reference genome assembly. Bacterial evolution is very different from either human or virus evolution, it's full of indels, missing genes, genes 'flying' around on plasmids and potential non-homologous recombination, risk of rearrangements etc ...

The challenge you face is automating the assemblies and then polishing those assemblies and there has been enormous progress recently, "closing the gaps" manually is a thing of the past - ancient history. There should be no under-estimation as the volume of recent work that has been dedicated particularly towards long read data. Even pre-assembly, pre-polishing the progress shouldn't be underestimated, albeit I suspect you're familiar with that. The old ideas of QUAST are now historic. The raised ante is noticeable.

Back to the original question consistently using a reference genome risks something going wrong that could cause an outlier artefact in the VCF and in my opinion its no longer cool for bacteria.

The rationale is clear, mis-assembly would cause an artefact in the tree and therefore in the epidemiological inference of the outbreak.

It's worth doing specific background reading because albeit I've deliberately not mentioned a given approach except for QUAST, the rate of progress has been significant. The world has changed.
> 1 votes

---
Tags: phylogenetics, variant-calling, read-mapping
---


title: Nextflow with Minimap2 not working

text:
# Question
Title: Nextflow with Minimap2 not working
I'm trying to create a Nextflow script that maps assemblies to reference gene sequences. This is what I have so far:

```

#!/usr/bin/env nexflow
nextflow.enable.dsl = 1

params.assmblreads = '/Users/rimo/Desktop/Assemblies/Assemblies_Flye/flye_assembly.fasta'
params.assmblref = '/Users/rimo/Desktop/Assemblies/gene_targets/GENE.fsa'

refgenome_file = file(params.assmblref)
fasta_files = Channel.fromPath(params.assmblreads)
                     .map { file -> tuple(file.baseName, file) }

process minimapping {

    input:
    file genome from refgenome_file
    set val(readID), file(readFile) from fasta_files

    output:
    set val(readID), file("${readID}.paf") into mapped_reads

    script:
    """
    ./minimap2 $readFile $genome > ${readID}.paf
    """

}

```

I'm running this on a macOS so I have to include the "./" before a command or it will give the error: command not found"

But everytime I run the script with the command: `./nextflow run trials.nf` it results in this error message:

```
(base) rimo@Ritas-MBP minimap2 % ./nextflow run ../attempt_out.nf
N E X T F L O W  ~  version 22.10.0
Launching `../attempt_out.nf` [nasty_baekeland] DSL1 - revision: 98d1deaebd
executor >  local (1)
[ae/47e3c0] process > minimapping (1) [  0%] 0 of 1
Error executing process > 'minimapping (1)'

Caused by:
executor >  local (1)
[ae/47e3c0] process > minimapping (1) [100%] 1 of 1, failed: 1 ✘
Error executing process > 'minimapping (1)'

Caused by:
  Process `minimapping (1)` terminated with an error exit status (1)

Command executed:

  ./minimap2 flye_assembly.fasta GENE.fsa > flye_assembly.paf

Command exit status:
  1

Command output:
  (empty)

Command error:
  .command.sh: line 2: ./minimap2: No such file or directory

Work dir:
  /Users/rimo/Desktop/AWS/minimap2/work/ae/47e3c0708f33b1949315fc83066b95

Tip: you can try to figure out what's wrong by changing to the process work dir and showing the script file named `.command.sh`

```

This is my first time ever writing a Nextflow script so I'm very new and confused to this new way of scripting. If someone could help me out that would be great!

# Answer
> 4 votes
> ```
> Command error:
>   .command.sh: line 2: ./minimap2: No such file or directory
>
> Work dir:
>   /Users/rimo/Desktop/AWS/minimap2/work/ae/47e3c0708f33b1949315fc83066b95
>
> ```

The problem is that your command is expecting the `minimap2` binary inside your working directory (i.e. the 'work dir', as defined above), but it could not be found. To avoid using relative (or absolute) paths, and to help keep your workflow portable, instead move your `minimap2` binary into a directory in your `$PATH`. Then, you should be able to call `minimap2` without a relative (or absolute) path to it. Here's one way using the new DSL2, which comes enabled by default with Nextflow 22.10.0:

```
params.assmblreads = '/Users/rimo/Desktop/Assemblies/Assemblies_Flye/flye_assembly.fasta'
params.assmblref = '/Users/rimo/Desktop/Assemblies/gene_targets/GENE.fsa'

process minimap2 {

    input:
    tuple val(readID), path(readFile)
    path genome

    output:
    tuple val(readID), path("${readID}.paf")

    """
    minimap2 "${readFile}" "${genome}" > "${readID}.paf"
    """
}

workflow {

    refgenome_file = file( params.assmblref )
    fasta_files = Channel.fromPath( params.assmblreads )
                         .map { file -> tuple( file.baseName, file ) }

    minimap2( fasta_files, refgenome_file )
}

```

To make your pipeline more portable, and to avoid manually installing and configuring your own tools, another way is to use Conda environments. Nextflow has built-in support for Conda and gives you access to all of the tools, including Minimap2, available in Bioconda. Install miniconda to get started, then add `conda.enabled = true` to your nextflow.config and the `conda` directive (e.g. `conda 'minimap2=2.24'`) to your process definition for example. In the same way, you could also use Docker (or Singularity) containers with the `container` directive. Note that Docker containers are also built by the Bioconda build system and these are uploaded to quay.io as part of the Biocontainers project (e.g. biocontainers/minimap2).

---
Tags: read-mapping, nextflow, errors, minimap2
---


# Question
Title: What does it mean when a gene and its transcript have opposite orientations in a GFF3 file?
I was working with a given GFF3 file, and I observed that some transcripts have orientation opposite to their transcripts. Here is a snippet:

```
chr1    .       gene    2189548 2194772 .       +       .       ID=c5eaf8c1-2bb2-566f-e154-6a9be797c11f;Name=FAAP20;gene_biotype=protein_coding;gene_id=34a96673-bfab-ebbf-b2b7-f9cb73163a4a;gene_name=FAAP20
chr1    .       transcript      2189548 2194772 .       -       .       ID=f47feeac-74ce-1543-53b6-afd9bbb8ab19;Name=rna-NM_182533.4;gene_biotype=protein_coding;gene_id=34a96673-bfab-ebbf-b2b7-f9cb73163a4a;gene_name=FAAP20;Parent=c5eaf8c1-2bb2-566f-e154-6a9be797c11f;transcript_biotype=protein_coding;transcript_name=rna-NM_182533.4;protein_id=NP_872339.3

```

How do I interpret this? I checked UCSC Browser and in each case it seems that the browser agrees with the transcript in each case, which seems like it should be canonical.

Moreover, I observed in this file that gene orientation is always `+`. Whereas a highly curated GFF3 of the human genome represents both directions:

```
# my file
$ awk '$3 == "gene" {print}' filtered_potential_genes.gff | cut -f7 | sort | uniq -c
   1600 +
# highly curated annotations
$ awk '$3 == "gene" {print}' /MANE.GRCh38.v0.95.select_refseq_genomic.gff | cut -f7 | sort|  uniq -c
   9393 +
   9191 -

```

According to the person who provided me the file, the existence of bidirectional transcription was indeed the motivation for this. However, it's obviously confusing to have it still be a definitive `+` rather than `.` or something to signify unstrandedness.

It seems likely that this is due to an artifact in my file. Are there other interpretations? Is the standard simply that gene orientation is meaningless and transcript orientation is important?

# Answer
> 2 votes
Looks like a poor annotation to me. I could understand a few scattered ncRNAs being arranged in the opposite direction from the gene, but all genes being oriented identically doesn't make sense.

---
Tags: sequence-annotation, gff3
---


# Question
Title: converting coordinates to sequences using bedtools getfasta - [segmentation fault] errors
I use **bedtools** to convert genomic coordinates (in the form of bed files, chrX x1:x2) to genomic sequence content in the form of fasta files ('ACGT...G').

The bedtools format for extracting sequence (*.fa file) from coordinates (*.bed file) is:

```
bedtools getfasta -fi  input_sequences.fa -bed  genomic_coordinates.bed -fo output_sequences.fa

```

Very frequently I get a "Segmentation fault (cored dumped)" during the run.

```
Warning: the index file is older than the FASTA file.
Segmentation fault (core dumped)

```

The strange thing is that the error will go away, if for example I add additional data (chromosomes) to the input fasta file (input\_sequences.fa), OR remove some chromosomes?? I am having trouble finding a logical path to explain this behavior so I can work around it.

There is no changing of fasta files while bedtools is working - in fact there is nothing in the primary directory other than the two files, the fasta file input\_sequence.fa and the bed file genomic\_coordinates.bed. There is 256GB of RAM memory on the machine, but that's not an issue, especially since the error goes away sometimes when additional chromosomes added to the fasta file.

Anyone experienced issues with bedtools getfasta erratic Segmentation faults? Known workarounds? I am using the current version of bedtools v2.30.0

What are some of the alternatives to "bedtools getfasta" for converting genomic coordinates to sequences?

Thanks...

# Answer
> 1 votes
The error about old index files is insightful to me. Bedtools presumably assumes that the input files are static / unchanging. If bedtools is automatically indexing the files before extracting sequences (when index files don't already exist), and the input files keep changing, that would be a good explanation for why problems are happening.

Do these input fasta files have bedtools and/or samtools index files in the same directory (usually `fileName.fa.fai`). Are you re-indexing the fasta files after changing them (e.g. via `samtools faidx fileName.fa`)?

I don't mean changing *during* the running of bedtools, I mean changing between different bedtools runs. If you create a file `input_sequences.fa`, then the first time it is subset, it will create an index file `input_sequences.fa.fai`. The next time bedtools is run to subset an `input_sequences.fa` within the same directory, it will check to see if an index file already exists; if one exists, it won't be regenerated. If the file `input_sequences.fa` has changed in the interim, then the index file will be incorrect, and the results of subsetting will be unpredictable.

This is why Bedtools reports the error: it can cause strange things to happen with output (intermittent segmentation faults, for example).

In fact, bedtools is explicit about this index generation in its output:

```
$ ls
sampleTagsNew.bed  sampleTagsNew.fa

$ bedtools getfasta -fi sampleTagsNew.fa -bed sampleTagsNew.bed
index file sampleTagsNew.fa.fai not found, generating...
>NewSampleTagX01:0-20
GTTGTCAAGATGCTACCGTT
>NewSampleTagX02:1-20
TTGTCAAGATGCTACCGTT

$ ls
sampleTagsNew.bed  sampleTagsNew.fa  sampleTagsNew.fa.fai

$ perl -i -pe 's/NewSampleTagX01/NewSampleTagX01_updated/' sampleTagsNew.fa

$ bedtools getfasta -fi sampleTagsNew.fa -bed sampleTagsNew.bed
Warning: the index file is older than the FASTA file.
>NewSampleTagX01:0-20
updatedGTTGTCAAGATG
>NewSampleTagX02:1-20
TagX02GTTGTCAAGATG

$ rm sampleTagsNew.fa.fai

$ bedtools getfasta -fi sampleTagsNew.fa -bed sampleTagsNew.bed
index file sampleTagsNew.fa.fai not found, generating...
WARNING. chromosome (NewSampleTagX01) was not found in the FASTA file. Skipping.
>NewSampleTagX02:1-20
TTGTCAAGATGCTACCGTT

```

In the example above, I carry out the following operations:

1. List the files in the directory
2. Run `bedtools getfasta` to subset a fasta file. Bedtools reports that an index file did not exist, so it is generating a new one.
3. List the files in the directory. Note that there is a new file: `sampleTagsNew.fa.fai`, which is the fasta file index.
4. Modify the sample tags fasta file \[to change the name of a sequence\]
5. Run `bedtools getfasta` again with an identical command line. It produces the same warning you have seen, and produces unexpected output, i.e. incorrect results.
6. Remove the fasta index file
7. Run `bedtools getfasta` again with an identical command line. This time the output is correct, stating that the first sequence cannot be found in the input fasta file.

---
Tags: bedtools, pybedtools
---


# Question
Title: How can I align two novel sequences to a reference genome and build a phylogenetic tree?
I have two new amino acid gene sequences that I want to align to a reference genome. Is it a good idea to have a combined phylogenetic tree? what will this type phylogenetic analysis called?

# Answer
> 2 votes
Edit: just to point out a tree needs 4 taxa (sequences), so rather than one reference genome, 3 reference genomes. Thus for Blast you might need to produce a database (using 3 reference genomes) before Blasting. There are other Blast approaches to obtaining a large meaningful sample size ... thats a separate question

---

It called a combined alignment, or concatenated genes.

The theory is simply that if two genes evolve at noticeably different rates they can't be combined into a single alignment and make a tree from them. The theory goes back a long but the modern test is in IQTREE.

You take your two alignments and test them for equal rates via IQTREE

```
iqtree2 -s example.phy -p example.nex --symtest-only

```

Testing for a combined alignment is also described in my past answer here: The previous question was can you combined an rRNA gene with a protein gene and the answer is ...

> Okay, you can't concatenate genes into a singular matrix particularly if these are rRNA genes versus protein genes, except in rare situations where the tests say 'thats okay' (below). This is simply because they evolve at very different rates rRNA genes have loci which evolve very slowly and are highly conserved.

---

The first part of the question was mapping the gene to a reference genome. Anything will do that ... local Blast works and using that is described in my answer to this question here:

If you loads of fasta seqs you'll need to make a database...

```
makeblastdb -in example700-1000.fa -parse_seqids -blastdb_version 5  -out laura700-1000.db -title "Laura 700 - 1000" -dbtype nucl

```

Then do this ...

```
blastn -query laura.fas -db laura700-1000.db -evalue 1e-6 -num_threads 4 -out blasthits.txt

```

However if its just reference vs. query then you can just use the above command.

---

The alignment is best done via a specialist alignment package Muscle Super5 and that is described in my answer here. This is done after your obtained the homologous regions in your reference genomes.

> HMM Super5 is a hidden Markov model (HMM) where gap penalties and the substitution matrix are the parameters at the centre of the model. That is the answer, thats why you can't change it - it does it itself and thats its purpose. Thus one algorithm (super5) is HMM for the gap-penalty and the other one I suspect isn't.

# Answer
> 0 votes
I would strongly recommend to find an outgroup or at least a fourth sequence. Building a phylogenetic tree with 3 sequences is doable but is tricky and error prone.

If you still want to do a 3-sequence analysis, you can align the two homologous proteins to the reference genome with exonerate, spaln or miniprot, extract the reference subsequence in the alignment, do a proper multiple-sequence alignment (MSA) with muscle or MAFFT, build a tree and put the root at the middle point of the longest branch (this assumes a molecular clock). After that, you can do typical phylogenetic analysis.

Whether the above works depends on the divergence of your sequences. If they are closely related and one sequence is a clear outgroup to the other two, the procedure could give you something, though having more sequences is still recommended.

EDIT in response to the comments: midpoint rooting assumes an approximate molecular clock across the whole tree. Severe violation of the assumption may lead to a wrong root. This approach is not very reliable due to the clock assumption.

---
Tags: python, phylogenetics, sequence-alignment, phylogeny
---


# Question
Title: How to visualize genome track of gene in specific cell-lines?
I'm trying to make a plot showing genome tracks of specific genes in specific cell-lines of RNA-seq and Chip-seq data.

It should look something like this I have recently seen this encode, but in this I didn't find any MDAMB231 cell-lines.

May I know how to make such plots.

I got the Counts data from CCLE portal.

```
             BT-549 MDAMB231
RP11-554I8.2    0   246
LINC00707      50   107
   RBM17    17036   10819
  PFKFB3    18098   7668
RP11-393K12.2   0   0

```

# Answer
You can use screen.encodeproject.org

It collected a lot of data from multiple cell lines or tissues. And you can export tracks to UCSC browser or other browsers.
> 1 votes

# Answer
*Answer from @devon-ryan, converted from comments:*

You'll need BAM or bigWig files for each of the tracks. If you can't find any you'll have to make them yourself from public fastq files. The fastq files can be obtained from ENA, SRA, or any of the other standard databases.
> 1 votes

# Answer
Like others have mentioned, these plots come from BAM, BED, and BigWig files. Those files contain the data that you can give to tools like the UCSC Genome Browser, IGV, HiGlass/Gosling, or other bioinformatics visualization software.

For the specific cell line you mention, you'll have to find the relevant RNA-seq and ChIP-seq data from some data repository, like CCLE, ENA, SRA, or GEO. You may be to put in some work converting available files to one of these compatible formats before plotting them.
> 1 votes

---
Tags: rna-seq, gene, genome, rna, encode
---


title: Error Registering Job Definition AWS

text:
# Question
Title: Error Registering Job Definition AWS
I am in the process of creating a bioinformatics pipeline and I am currently trying to setup the AWS infrastructure for which it will run.

The script I am running to setup the AWS batch job looks like the following:

```
#! /bin/bash

# Setup the AWS Batch architecture
# Copyright 2022, SolareaBio

set -e

# ***** BEGIN custom lines

container="BGC_bacteria"
vCPUs_req=12
memory_MiB_req=40000
# ASK ABOUT CUSTOM CPU AND MEMORY SETTINGS

# ***** END custom lines

container_lc=$(echo $container | tr [A-Z] [a-z])

subnet_useast1a="subnet-42371d0a"
subnet_useast1b="subnet-b03c54ea"
subnet_useast1c="subnet-00ee1364"

usage() {
  cat << EOF
Usage $0 [-h]

Builds the $container docker container, pushes it to the AWS ECR, and sets up the $container AWS Batch environment
EOF
}

while getopts "h" option; do
  case $option in
    h)   usage;exit;;
    [?]) usage; exit;;
  esac
done

script_dir=$(dirname "$0")

repo="990850896279.dkr.ecr.us-east-1.amazonaws.com"

echo "Logging into ECR"
aws ecr get-login-password --region us-east-1 | docker login --username AWS --password-stdin ${repo}

image_and_tag="${container_lc}:latest"
echo "Building the docker container: $image_and_tag"

# If you get errors from the docker build command, try adding --no-cache
docker build --no-cache -f $script_dir/${container}.dockerfile -t $image_and_tag $script_dir

echo "Tagging Image"
docker tag ${image_and_tag} ${repo}/${image_and_tag}

# Create repository in AWS ECR
repo_container_info=$(aws ecr describe-repositories --repository-names ${container_lc})
if [[ "$repo_container_info" == "" ]]; then
  aws ecr create-repository --repository-name ${container_lc}
fi

echo "Pushing Image up to AWS"
docker push ${repo}/${image_and_tag}

# Create compute environment
EC2_type=ON_DEMAND
ce_info=$(aws batch describe-compute-environments --compute-environments ${container}_${EC2_type}_CE)
ce_status=$(echo $ce_info | jq '.computeEnvironments[0] | .status // ""' | tr -d '"')
if [[ $ce_status != "VALID" ]]; then
  compute_resources='{"type":"EC2","allocationStrategy":"BEST_FIT_PROGRESSIVE","minvCpus":0,"maxvCpus":10000,"tags":{"Name":"'${container}_Batch'"},"instanceTypes":["c5a","m5a","c5","m5"],"subnets":["'$subnet_useast1a'","'$subnet_useast1b'","'$subnet_useast1c'"],"securityGroupIds":["sg-ff2fd48f","sg-65ff5315"],"ec2KeyPair":"SBTools","instanceRole":"ecsInstanceRole"}'
  aws batch create-compute-environment \
    --compute-environment-name ${container}_${EC2_type}_CE \
    --type MANAGED \
    --state ENABLED \
    --compute-resources $compute_resources \
    --service-role AWSBatchServiceRole

  printf "Waiting for Compute Environment to be created"
  while [[ $ce_status != "VALID" ]]; do
sleep 5
printf "."
ce_info=$(aws batch describe-compute-environments --compute-environments ${container}_${EC2_type}_CE)
    ce_status=$(echo $ce_info | jq '.computeEnvironments[0] | .status // ""' | tr -d '"')
  done
  echo ""
fi

# Create the job queue
jq_info=$(aws batch describe-job-queues --job-queues ${container}_${EC2_type}_JQ)
jq_status=$(echo $jq_info | jq '.jobQueues[0] | .status // ""' | tr -d '"')
if [[ $jq_status != "VALID" ]]; then
  aws batch create-job-queue \
    --job-queue-name ${container}_${EC2_type}_JQ \
    --state ENABLED \
    --priority 5 \
    --compute-environment-order '{"order":1,"computeEnvironment":"'${container}_${EC2_type}_CE'"}'
fi

# Register job definition
job_definition="${container}_JD"
jd_info=$(aws batch describe-job-definitions --status ACTIVE --job-definition-name $job_definition)
jd_status=$(echo $jd_info | jq '.jobDefinitions[0] | .status // ""' | tr -d '"')
if [[ $jd_status != "ACTIVE" ]]; then
  aws batch register-job-definition \
    --job-definition-name $job_definition \
    --type container \
    --container-properties '{"privileged":true,"image":"'$repo/$image_and_tag'","vcpus":'${vCPUs_req}',"memory":'${memory_MiB_req}',"mountPoints":[{"sourceVolume":"SB-EFS","containerPath":"/EFS"}],"volumes":[{"host":{"sourcePath":"/"},"name":"SB-EFS"}],"command":["Ref::container","Ref::run_id","Ref::input_file"]}' \
    --retry-strategy '{"attempts": 10,"evaluateOnExit":[{"onStatusReason" :"Host EC2*","action": "RETRY"},{"onReason" : "*","action": "EXIT"}]}'
fi

# Make sure the EFS has an access point set up
efs_id=$(aws efs describe-file-systems | jq -r '.FileSystems[] | select (.Name=="SB-EFS") | .FileSystemId')
if [[ "$efs_id" == "" ]]; then
  echo "ERROR: SB-EFS does not exist"
  exit
fi
ap_id=$(aws efs describe-access-points | jq -r '.AccessPoints[] | select (.Name=="EFS_AP") | .AccessPointId')
if [[ "$ap_id" == "" ]]; then
  # 1001:1001 is user sb in group sb
  ap_id=$(aws efs create-access-point --file-system-id "$efs_id" --posix-user '{"Uid":1001,"Gid":1001}' --root-directory '{"Path":"/"}' --tags '[{"Key":"Name","Value":"EFS_AP"}]' | jq -r '.AccessPointId')
fi

```

Now when I run it, everything works as it should until getting to the --retry-strategy command line in the "register job definition" section of the script. The following error is subsequently thrown:

```
Parameter validation failed:
Unknown parameter in retryStrategy: "evaluateOnExit", must be one of: attempts

```

Not sure what the issue is with the "attempts" parameter that is set. Any potential solutions or guidance to solving this error is much appreciated. Thanks!

# Answer
> 0 votes
You didn't say what version of the AWS CLI you're using, but I suspect it's version 1. It looks like the `evaluateOnExit` feature was added in version 2 but was not backported. Installing AWS CLI version 2 resolves the issue:

https://docs.aws.amazon.com/cli/latest/userguide/getting-started-install.html

Note that Nextflow also provides built-in support for AWS Batch. You'd still need to configure some AWS infrastructure using Nextflow, but it should be a lot less painful than writing/maintaining your own script(s). It's definitely worth checking out if you haven't already.

---
Tags: bash, shell, aws
---


# Question
Title: Modification time does not trigger a snakemake rule
I'm trying to use snakemake to replace some make workflows, but I'm finding that snakemake does not use the file modification time to trigger a rule.

Here is my snakefile

```
rule sample:
    input:
        "a"
    output:
        "b"
    shell:
        "touch {output}"

```

If I create `a` then run snakemake, `b` is created as expected:

```
$ touch a
$ snakemake -c1 b
Building DAG of jobs...
Using shell: /usr/bin/bash
Provided cores: 1 (use --cores to define parallelism)
Rules claiming more threads will be scaled down.
Job stats:
job      count    min threads    max threads
-----  -------  -------------  -------------
sample       1              1              1
total        1              1              1

Select jobs to execute...

[Sat Oct 29 14:36:43 2022]
rule sample:
    input: a
    output: b
    jobid: 0
    reason: Missing output files: b
    resources: tmpdir=/tmp

[Sat Oct 29 14:36:43 2022]
Finished job 0.
1 of 1 steps (100%) done
Complete log: .snakemake/log/2022-10-29T143643.334050.snakemake.log
$ ls -l a b
-rw-r--r-- 1 ... ... 0 Oct 29 14:36 a
-rw-r--r-- 1 ... ... 0 Oct 29 14:36 b

```

Now, if I touch `a` then ask snakemake to update `b`, nothing happens:

```
$ touch a
$ snakemake -c1 b
Building DAG of jobs...
Nothing to be done (all requested files are present and up to date).
Complete log: .snakemake/log/2022-10-29T144005.589943.snakemake.log
$ ls -l a b
-rw-r--r-- 1 ... ... 0 Oct 29 14:40 a
-rw-r--r-- 1 ... ... 0 Oct 29 14:36 b

```

Modifying the contents of `a` triggers an update as expected:

```
$ echo new > a
$ snakemake -c1 b
Building DAG of jobs...
Using shell: /usr/bin/bash
Provided cores: 1 (use --cores to define parallelism)
Rules claiming more threads will be scaled down.
Job stats:
job      count    min threads    max threads
-----  -------  -------------  -------------
sample       1              1              1
total        1              1              1

Select jobs to execute...

[Sat Oct 29 14:42:13 2022]
rule sample:
    input: a
    output: b
    jobid: 0
    reason: Updated input files: a
    resources: tmpdir=/tmp

[Sat Oct 29 14:42:13 2022]
Finished job 0.
1 of 1 steps (100%) done
Complete log: .snakemake/log/2022-10-29T144213.161985.snakemake.log
$ ls -l a b
-rw-r--r-- 1 ... ... 4 Oct 29 14:42 a
-rw-r--r-- 1 ... ... 0 Oct 29 14:42 b

```

I installed snakemake-minimal (version 7.17.1) from bioconda, am using python 3.9.13 from conda-forge, on a RHEL 7 system.

```
$ snakemake --version
7.17.1
$ python --version
Python 3.9.13
$ cat /etc/redhat-release
Red Hat Enterprise Linux Server release 7.8 (Maipo)

```

I'm clearly misunderstanding something very basic about snakemake. Any suggestions on how to debug this, descriptions of how to have snakemake pay attention to file modification times, or just an explanation of why snakemake behaves differently from make would be much appreciated!

thanks in advance,

# Answer
This is a Snakemake feature: For small files (less than 100000 bytes), checksums are compared to determine if a job needs to run if the input file is newer than output file:

https://github.com/snakemake/snakemake/pull/1568

With a size of 99999 bytes:

```
$ dd if=/dev/zero of=a bs=99999 count=1
1+0 records in
1+0 records out
99999 bytes (100 kB, 98 KiB) copied, 0.000502722 s, 199 MB/s
$ snakemake -c1 b
Building DAG of jobs...
Using shell: /usr/bin/bash
Provided cores: 1 (use --cores to define parallelism)
Rules claiming more threads will be scaled down.
Job stats:
job       count    min threads    max threads
------  -------  -------------  -------------
sample        1              1              1
total         1              1              1

Select jobs to execute...

[Mon Oct 31 14:38:00 2022]
rule sample:
    input: a
    output: b
    jobid: 0
    reason: Missing output files: b
    resources: tmpdir=/tmp

[Mon Oct 31 14:38:00 2022]
Finished job 0.
1 of 1 steps (100%) done
Complete log: .snakemake/log/2022-10-31T143800.486966.snakemake.log
$ dd if=/dev/zero of=a bs=99999 count=1
1+0 records in
1+0 records out
99999 bytes (100 kB, 98 KiB) copied, 0.000180857 s, 553 MB/s
$ snakemake -c1 b
Building DAG of jobs...
Nothing to be done (all requested files are present and up to date).
Complete log: .snakemake/log/2022-10-31T143813.285988.snakemake.log

```

With a size of 100000 bytes:

```
$ dd if=/dev/zero of=a bs=100000 count=1
1+0 records in
1+0 records out
100000 bytes (100 kB, 98 KiB) copied, 0.000173008 s, 578 MB/s
$ snakemake -c1 b
Building DAG of jobs...
Using shell: /usr/bin/bash
Provided cores: 1 (use --cores to define parallelism)
Rules claiming more threads will be scaled down.
Job stats:
job       count    min threads    max threads
------  -------  -------------  -------------
sample        1              1              1
total         1              1              1

Select jobs to execute...

[Mon Oct 31 14:41:14 2022]
rule sample:
    input: a
    output: b
    jobid: 0
    reason: Missing output files: b
    resources: tmpdir=/tmp

[Mon Oct 31 14:41:14 2022]
Finished job 0.
1 of 1 steps (100%) done
Complete log: .snakemake/log/2022-10-31T144114.211255.snakemake.log
$ dd if=/dev/zero of=a bs=100000 count=1
1+0 records in
1+0 records out
100000 bytes (100 kB, 98 KiB) copied, 5.35361 s, 18.7 kB/s
$ snakemake -c1 b
Building DAG of jobs...
Using shell: /usr/bin/bash
Provided cores: 1 (use --cores to define parallelism)
Rules claiming more threads will be scaled down.
Job stats:
job       count    min threads    max threads
------  -------  -------------  -------------
sample        1              1              1
total         1              1              1

Select jobs to execute...

[Mon Oct 31 14:41:26 2022]
rule sample:
    input: a
    output: b
    jobid: 0
    reason: Updated input files: a
    resources: tmpdir=/tmp

[Mon Oct 31 14:41:26 2022]
Finished job 0.
1 of 1 steps (100%) done
Complete log: .snakemake/log/2022-10-31T144125.949828.snakemake.log

```
> 5 votes

---
Tags: snakemake
---


# Question
Title: How to adjust #columns shown when PyRange is printed?
In Pandas, we have `pandas.options.display` to adjust how a DataFrame should be displayed when printed.

I wonder if `PyRanges` has a similar set of options?

For example, if I have a PyRange object like this:

```
import pyranges as pr
data = pr.data.ucsc_bed()

```

I can do `print(data)` and get:

```
+--------------+-----------+-----------+------------+------------+-------+
| Chromosome   | Start     | End       | Feature    | gene_id    | +4    |
| (category)   | (int32)   | (int32)   | (object)   | (object)   | ...   |
|--------------+-----------+-----------+------------+------------+-------|
| chr1         | 12776117  | 12788726  | gene       | AADACL3    | ...   |
| chr1         | 169075927 | 169101957 | gene       | ATP1B1     | ...   |
| chr1         | 6845383   | 7829766   | gene       | CAMTA1     | ...   |
| chr1         | 20915589  | 20945396  | gene       | CDA        | ...   |
| ...          | ...       | ...       | ...        | ...        | ...   |
| chrX         | 152661096 | 152663330 | exon       | PNMA6E     | ...   |
| chrX         | 152661096 | 152666808 | transcript | PNMA6E     | ...   |
| chrX         | 152664164 | 152664378 | exon       | PNMA6E     | ...   |
| chrX         | 152666701 | 152666808 | exon       | PNMA6E     | ...   |
+--------------+-----------+-----------+------------+------------+-------+
Stranded PyRanges object has 5,519 rows and 9 columns from 30 chromosomes.
For printing, the PyRanges was sorted on Chromosome and Strand.
4 hidden columns: transcript_id, Strand, exon_number, transcript_name

```

How can I increase the number of columns shown?

Also appreciated if you refer me to a documentation that explains other related display options.

**Update:** For future reference, @terdon remark encouraged me to check if such an adjustment is even considered in the `PyRanges`. I found the line that is responsible for this representation, and indeed the `COLUMN` environment variable is the only defining parameter. In conclusion, the adjustment that I am looking for is simply not implemented in `PyRanges`.

# Answer
When trying this on my machine, if I resize the terminal window before running `print(data)`, I get the full output:

You can also force this, even if your terminal is actually smaller, by setting the `COLUMNS` variable to something larger than the actual size when layching `python` or `ipython`. This is a special bash variable that holds the current size of the terminal. From `man bash`:

```
   COLUMNS
          Used  by  the select compound command to determine the terminal
          width when printing selection lists.  Automatically set if  the
          checkwinsize  option is enabled or in an interactive shell upon
          receipt of a SIGWINCH.

```

So if you set it manually, you can get the entire data frame printed, BUT you don't really want to do that since it won't scroll, it just wraps the text. For example, try:

```
COLUMNS=300 python

```

That will set the size to 300, which results in this:
> 1 votes

---
Tags: pyranges
---


# Question
Title: How to find closely related genes for a specific gene from WGCNA modules?
I have used WGCNA by integrating several datasets which are processed in a similar way. Identified 17 modules, followed by performed pathway analysis with the genes present in those modules.

Among all the protein-coding genes,ncRNAs and other genes, I'm particularly interested in a specific lncRNA `MALAT1`.

I found that `MALAT1` is co-expressed with more than 800 genes in a module named `blue`. I have few questions, can anyone please help me in this?

Questions:

1. I would like to check which are the closest co-expressed genes to `MALAT1`. How can I know that?
2. Once I get the closest genes to `MALAT1`, I would like to make a network with those genes and MALAT1. How to do this network?

# Answer
> 3 votes
The simplest way is to re-calculate (or load, if saved) the TOM matrix and find the genes with the highest TOM to MALAT1. If you have enough RAM to carry out WGCNA in a single block, you can use TOMsimilarityFromExpr, make sure you give it the same arguments as the ones you used in network construction, and select the appropriate column from the result. Then order the column and retain top genes, as many as you think is reasonable. You can also use adjacency instead of TOM, though TOM is a more natural quantity to use if you used in your network construction. Once you figure out which genes you want to retain, you can subset the TOM matrix (both rows and columns) to those genes and you have the sub-network among these genes.

---
Tags: r, rna-seq, networks, wgcna
---