title: Singularity and Illumina's Nirvana

text:
# Question
Title: Singularity and Illumina's Nirvana
I am trying to using Illumina's Nirvana however I can't seem to get it to work despite having the right Singularity path.

1. I pulled Nirvana down using

```
singularity pull docker://annotation/nirvana:3.14

```
2. I found the path to the .dll using the following command:

```
singularity exec nirvana_3.14.sif /bin/bash -c "cd /opt/nirvana/ && ls"

```
3. Tried to run the downloader using:

```
## Define path to singularity image:
sif_container="/path/to/nirvana_3.14.sif"

## Path in the singularity file:
sif_path="/opt/nirvana/Downloader.dll"

## Data to be placed in
Data="path/to/Nirvana/GRCh38/"

#_____________________________________________________________________________________________________________________________________

singularity exec "$sif_container" dotnet "$sif_path" --ga GRCh38 --out "$Data"

```

yet it keeps saying `ERROR: The top-level output directory directory was not specified. Please use the --out parameter.`. I have checked the paths and they are correct and I have specified the data path yet it does not seem to recognize it? I have even moved around the container and the path variables to no avail?

# Answer
> 3 votes
From your previous questions, it is apparent that your HPC has an odd singularity setup. You could ask you sysadmin about how to navigate it.

## EDIT

After working with OP on chat (and based on previous experience with them), it looks like their sysadmin has `mount hostfs` enabled in the `singularity.conf` file, which messes with every container every time. I guess it's easier than letting people bind their directories of choice but anyway, the way to circumvent that is to use `--no-mount hostfs`

@Indira, always run `singularity [exec|shell|run] --no-mount hostfs` instead of plain `singularity [exec|shell|run]`. All containers will then run well and if they error out, it's probably because either the container internally encountered an error or it was unable to access an external resource. What you're running into now is a mix that obscures the real problems and creates a no-win scenario for you.

---
Tags: illumina, singularity
---


# Question
Title: How to calculate the variant allelic fraction (VAF) from VEP output?
*This question was also asked on Biostars*

How can I calculate the variant allelic fraction from the VEP file for multiple human samples?

I used the ensembl-vep and it is human data.

# Answer
> 2 votes
If you just have a VEP output (i.e. VCF) without raw genotypes then you cannot calculate the allelic fraction.

---
Tags: python, biopython, vep
---


title: Reading same BAM file twice with htslib

text:
# Question
Title: Reading same BAM file twice with htslib
I would like to read iteratively the same BAM file twice using `htslib`.

My option is to use `hts_open`, `sam_hdr_read`, then `sam_read1` as much as needed to read the file until EOF.

Is there a possibility to "rewind", and use `sam_read1` again from the beginning of the file, without closing the file, destroying the object, and opening it again?

---

**EDIT**

My aim is group the reads by the barcode. The first pass would allocate the specific sizes first (the number of reads by barcode), then store the reads.

I know I could sort the reads first (using `samtools sort -t BX`), but I would like to skip this step using pre-allocation.

# Answer
> 1 votes
The best solution is to reopen the file. For a coordinate sorted file you have the option to jump back to the start of the first chromosome with index but that is more complex, less flexible and won't be faster.

---
Tags: htslib
---


title: Residue numbering in biopython

text:
# Question
Title: Residue numbering in biopython
Is there an out-of-the-box function in the `Bio.SeqIO` module (in `biopython`) that can print a protein sequence (in command line) in multiple lines along with the numbering of the residues (either left and right or in-between the lines)?

In particular I have a `.fasta` file which I can read into a `SeqRecord` object like so:

```
from Bio import SeqIO
record = SeqIO.read("protein.fasta", "fasta")

```

I would like to then print the `Seq` class with the numbering, not just as a string as shown below:

```
print(record.seq)

```

I have already tried:

```
SeqIO.write(record, sys.stdout, 'gb')

```

But it is not working for me. The error was `ValueError: missing molecule_type in annotations`. I googled it and it appears to be related to newer biopython versions (see https://github.com/linsalrob/PhiSpy/issues/43).

This on-line tool seems to do what I want, but one needs to provide it with an alphabet. How do I do that with code? Does a `SeqRecord` object have an alphabet?

# Answer
> 3 votes
@M\__ is right - there's a `molecule_type` error when you try to write and this needs to be fixed. Here's my end to end approach:

The input file:

```
cat xyz.fa
>NP_000537.3 cellular tumor antigen p53 isoform a [Homo sapiens]
MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGPDEAPRMPEAA
PPVAPAPAAPTPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGFLHSGTAKSVTCTYSPALNKMFCQLAKT
CPVQLWVDSTPPPGTRVRAMAIYKQSQHMTEVVRRCPHHERCSDSDGLAPPQHLIRVEGNLRVEYLDDRN
TFRHSVVVPYEPPEVGSDCTTIHYNYMCNSSCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGR
DRRTEEENLRKKGEPHHELPPGSTKRALPNNTSSSPQPKKKPLDGEYFTLQIRGRERFEMFRELNEALEL
KDAQAGKEPGGSRAHSSHLKSKKGQSTSRHKKLMFKTEGPDSD

```

My (bad) code:

```
from Bio import SeqIO

with open("xyz.fa") as handle:
    for record in SeqIO.parse(handle, "fasta"):
        record.annotations = { "molecule_type" : "protein" }
SeqIO.write(record, "xyz.gb", "gb") # moved two indentations to the left

```

The output:

```
cat xyz.gb

LOCUS       NP_000537.3              393 aa                     UNK 01-JAN-1980
DEFINITION  NP_000537.3 cellular tumor antigen p53 isoform a [Homo sapiens].
ACCESSION   NP_000537
VERSION     NP_000537.3
KEYWORDS    .
SOURCE      .
  ORGANISM  .
            .
FEATURES             Location/Qualifiers
ORIGIN
        1 meepqsdpsv epplsqetfs dlwkllpenn vlsplpsqam ddlmlspddi eqwftedpgp
       61 deaprmpeaa ppvapapaap tpaapapaps wplsssvpsq ktyqgsygfr lgflhsgtak
      121 svtctyspal nkmfcqlakt cpvqlwvdst pppgtrvram aiykqsqhmt evvrrcphhe
      181 rcsdsdglap pqhlirvegn lrveylddrn tfrhsvvvpy eppevgsdct tihynymcns
      241 scmggmnrrp iltiitleds sgnllgrnsf evrvcacpgr drrteeenlr kkgephhelp
      301 pgstkralpn ntssspqpkk kpldgeyftl qirgrerfem frelnealel kdaqagkepg
      361 gsrahsshlk skkgqstsrh kklmfktegp dsd
//

```

I used the code from https://github.com/biopython/biopython/issues/3156 to fix the `ValueError: missing molecule_type in annotations` that I faced when I'd skipped the `record.annotations = {"molecule_type" : "protein" }` line.

There is definitely a better way to write the python code instead of using a loop or at least writing to a static file from within a loop, but that's not central to the problem here so I'm choosing to ignore it. I welcome any improvements to the code.

# Answer
> 0 votes
Given the criteria required the closest solution is:

```
from Bio import SeqIO

seq_obj = ... # Biopython seq_obj
with open("mysequence.gb", "w") as fout:
    SeqIO.write(seq_obj, fout, "gb")

```

If this fails I would substitute `gb` for `fasta` and try again. If that fails, i.e. no fasta, there's a problem with the sequence object. If this succeeds but the `gb` fails, just supply a reproducible example. It will work.

I've used `gb` to successfully write APIs against NCBI Genbank and trap the meta-data. Genbank format is really cool.

Oh wait, sorry the answer is ...

```
seq_obj for records:
    print(seq_obj.format("gb"))

```

---

From the comments "molecule type"...

Which are you using

```
for record in SeqIO.parse(filename, "fasta"):

```

or

```
for record in SeqIO.read(filename, "fasta"):

```

Which ever is used I'd to the other method.

Just to comment on @RamRS solution, thats cool I didn't know.

```
for record in SeqIO.parse(filename, "fasta"):
    record.annotations = {"molecule_type" : "protein" }

```

What is happening is the instance is being manually set. However, Biopython should detected the data type automatically. My suspicion is that `read` or `parse` one comprises and automated data type setter.

In any case, RamRS's solution is useful and will work. In past coding I would have read the object via Genbank.

The old `alphabet` method is deprecated and shouldn't be attempted.

---

For future reference RamRS's code is described here

---
Tags: sequence, biopython, command-line
---


title: Gviz Coverage Plots

text:
# Question
Title: Gviz Coverage Plots
*This question has also been asked on Biostars*

I have two bam files from single-cell RNA sequencing mapped to the reference genome using CellRanger, I can view them in IGV and I have a particular region where the pattern of reads mapped to the reference genome are different between the two bam files but when I try the following code to plot trackplots with coverage, the resulting plot is not consistent with what I see in IGV. I am not very familiar with this, aren't I supposed to normalize the data as well? please help me.

```
txdb <-TxDb.Hsapiens.UCSC.hg38.knownGene
txTr <- GeneRegionTrack(txdb, chromosome = "chr1",
                        start = 153561336,  end = 153561564)
dTrack <- DataTrack(range = bam.file, genome = "hg38", type = "l",
                     name = "Coverage", window = -1,
                     chromosome = "chr1")

dTrack0 <- DataTrack(range = bam.file.0, genome = "hg38", type = "l",
                     name = "Coverage", window = -1,
                     chromosome = "chr1", col="red")

plotTracks(c(dTrack,dTrack0, txTr), from = 153561357, to = 153561585)

```

# Answer
> 0 votes
I would recommend reading the following study by Zhao *et al* (2021).

The traditional simplest measure is CPM otherwise known as TPM. This is counts/transcripts per million. FPKM (fragments per kilobase of transcript per million fragments mapped) is used for metagenomics. I don't know its application in your context ... I don't know your context.

The discussions for the pros- and cons- of FPKM and TPM/CPM are discussed therein.

The authors in the study favour normalisation - *in their context* \- not in general. The authors use the coefficient of variation, and intraclass correlation coefficient. I would see this as a form of standardisation. It's worth keeping an open mind because its a pivotal step in the calculation but all too easily overlooked area of investigation.

From RamRS's comments that may have changed already. In general statistics standardisation is a pivotal step, usually followed by the removal of outliers.

---
Tags: r, scrnaseq, sam, cellranger, gviz
---


# Question
Title: How can I detect the Apo B gene mutations in my raw dna data?
Familial defective apolipoprotein B-100 (FDB) is an autosomal dominant disorder associated with increased plasma cholesterol concentration and may thus increase the risk of premature coronary heart disease. To date, three apo B mutations, apo B R3500Q, apo B R3531C and apo B R3500W that cause defective apo B, have been reported.

I have my raw DNA data from 23andMe (where each line corresponds to a single SNP). How can I search it for those related Apo B gene mutations?

A single line in my DNA data looks like this:

```
rs12567639 1 756268 AA

```

I have a programming background, but none in the bioinformatics area.

# Answer
> 3 votes
I think you mean these. The amino acid numbering has changed, but that happens as we discover more transcript variants.

http://useast.ensembl.org/Homo\_sapiens/Variation/Mappings?db=core;r=2:21005788-21006788;v=rs5742904;vdb=variation;vf=184237596

http://useast.ensembl.org/Homo\_sapiens/Variation/Mappings?db=core;r=2:21005696-21006696;v=rs12713559;vdb=variation;vf=186851457

http://useast.ensembl.org/Homo\_sapiens/Variation/Explore?db=core;g=ENSG00000084674;r=2:21001429-21044073;v=rs144467873;vdb=variation;vf=251617275

# Answer
> 2 votes
If you are interested in ApoB gene variants you can:

https://www.ensembl.org/Homo\_sapiens/Variation/Explore?db=core;g=ENSG00000084674;r=2:21001429-21044073;v=rs5742904;vdb=variation;vf=184237596

But I assume all of this is moot since 23andMe should give you some summary of the results and not just the raw file, right?

---
Tags: snp
---


# Question
Title: Identification of unique chains in proteins by sequence similarity
I'm trying to find unique chains or proteins within their PDB files. Many proteins have multiple chains, but very often they are identical (say, the PDB file consists of 4 homodimers, for example, so there is one unique chain represented 8 times total). I am processing certain structural features of the chains, so processing two identical chains is redundant for my pruposes.

I am using Biopython and the way I'm currently identifying unique chains is looking through each PDB file, checking the sequences of the chains and comparing them to one another using alignments (specifically `pairwise2.align.globalxx(seq1, seq2)`. I consider a chain non-unique if it shares \>90% sequence similarity to another chain. (Most chains have 100+ residues.)

My question is, is this threshold appropriate? Should I go lower, maybe 80%, or the opposite, higher? For now, it seems to work fine, but I'll be processing a lot of PDB files and want to be on the safer side. If I accidentally process identical chains, that's not much of an issue, but it adds unnecessary time to the processing. Is there maybe a **source/article** to help with what sequence similarity I should be using? I've found multiple articles, but they usually talk about homology (i.e., common ancestry of proteins), but I'm only doing this within the protein itself and its subunits (chains), so this is of no use for me.

# Answer
> 3 votes
Came across this looking for a solution myself last night but now I have created my own I will share it.

Rather then using `Bio`python (which is a Python package I use a lot), it is more computationally efficient to parse the `.fasta` files from the protein data bank (RCSB in this case) which as far as I can see, contain this information. For example, this is 2P16's `.fasta` file (with 2 unique chains):

```
>2P16_1|Chain A|Coagulation factor X (EC 3.4.21.6) (Stuart factor) (Stuart-Prower factor)|Homo sapiens (9606)
IVGGQECKDGECPWQALLINEENEGFCGGTILSEFYILTAAHCLYQAKRFKVRVGDRNTEQEEGGEAVHEVEVVIKHNRFTKETYDFDIAVLRLKTPITFRMNVAPACLPERDWAESTLMTQKTGIVSGFGRTHEKGRQSTRLKMLEVPYVDRNSCKLSSSFIITQNMFCAGYDTKQEDACQGDSGGPHVTRFKDTYFVTGIVSWGEGCARKGKYGIYTKVTAFLKWIDRSMKT
>2P16_2|Chain B[auth L]|Coagulation factor X (EC 3.4.21.6) (Stuart factor) (Stuart-Prower factor)|Homo sapiens (9606)
KLCSLDNGDCDQFCHEEQNSVVCSCARGYTLADNGKACIPTGPYPCGKQTLE

```

and 6WON's (with 8 of the same chains):

```
>6WON_1|Chains A, B, C, D, E, F, G, H|YohF|Salmonella typhimurium (strain LT2 / SGSC1412 / ATCC 700720) (99287)
QGMTKVAIVTASDSGIGKACALLLAQNGFDIGITWHSDERGAQETAKKAAQFGVRAETIHLDLSQLPEGAQAIEHLIQRLGRVDVLVNNAGAMTKSAFIDMPFTQWRQIFTVDVDGAFLCAQIAARHMIKQGEGGRIINITSVHEHTPLPQASAYTAAKHALGGLTKSMALELIEHHILVNAVAPGAIATPMNDMDDSDIEPGSEPSIPIARPGSTHEIASLVAWLCSEGASYTTGQSLIVDGGFMLANPQFNAK

```

So far this always seems to be the case and parsing the `.fasta` is a lot easier then using BLAST or sequence similarity such as you are or were.

For the implementation, whereas the `.pdb`, `.cif` files etc use the same download URL (`https://files.rcsb.org/download/`), the `.fasta` files use this link: `https://www.rcsb.org/fasta/entry/`. Parsing the downloaded `.fasta` files and only reading the lines that begin with `>`, stripping out unnecessary chains (in my specific case I just want the first unique chains in the PDB so `A, B` being the same and `C, D` being the same would return me a list with \[`A, C`\], you can modify this script to your own ends) and running the Python code gives this output:

```
First unique chains in 1d2z (2 unique chains): ['B', 'A']
First unique chains in 2p16 (2 unique chains): ['A', 'B']
First unique chains in 6won (1 unique chain):  ['A']

```

There was an issue with some entries having `Chain` or `Chains` so both had to be considered. Additionally, there are commas that can be removed from all but the last chain ID but for my purposes I do not need to strip it, you can put something in `_fasta_helper()` to this end before `return`. This is the Python code and hopefully it helps someone out in the future:

```
import re
from typing import Union
import urllib.request
from urllib.error import (HTTPError,
                          URLError,
                          ContentTooShortError as ContentError)

def rcsb_downloader(pdb: str, ext: str = None, down_path: str = None) -> str:
    """Downloads the PDB from the protein databank RCSB.

    Parameters
    ----------
    pdb : str
        Protein databank accession 4 letter code i.e. '6won'
    ext : str, default = None
        Extension of the file to download i.e. '.fasta', '.cif' etc. Note: if
        None argument then ext = '.pdb'

    Returns
    -------
    str | bool
        Returns str path if download successful, False if failure
    """

    if ext is None:
        ext = '.pdb'
    if down_path is None:
        download_path = ('/your/path/to/dir/' + pdb + ext)
    else:
        download_path = down_path + pdb + ext

    rcsb_path = 'https://files.rcsb.org/download/' + pdb + ext

    if ext == '.fasta': # fasta has a different web search then others
        rcsb_path = 'https://www.rcsb.org/fasta/entry/' + pdb

    try:
        urllib.request.urlretrieve(rcsb_path, download_path)
        return download_path
    # Try again in case there was a temporary issue with the first attempt
    except (ContentError, HTTPError, URLError):
        try:
            urllib.request.urlretrieve(rcsb_path, download_path)
            return download_path
        except (ContentError, HTTPError, URLError):
            print(f"PDB failed to download {pdb}")
            return False

    def _fasta_helper(self, line: str, single_unique_chain: bool = True) -> Union[str, list]:
        """Separate the first chain ID from the .fasta comment line."""

        # [1] some of the chains had [auth C] etc after so this is removed
        # [2] was originally a list comprehension that filtered based on the word
        # chain, however some of the fasta entries had "Chain" or "Chains" in their
        # description which broke this. The first arrow is removed, split at |, the
        # second entry (1st in Python) entry contains the chain information and the
        # first entry removed (which so far has always been chain, anticipate their
        # being edge cases with older PDB entries). This appears to hold up.

        clean_chain = [
            x.replace(",", "")
            for x in list(
                re.sub(r"\[.*?\]", "", line).lstrip(">").split("|")[1].split(" ")[1:]  # [1]  # [2]
            )
        ]

        if not single_unique_chain:
            return clean_chain
        return clean_chain[0]

def fasta_parse(path_: str) -> list:
    """Parse .fasta files and analyse lines beginning with '>'."""

    with open(path_, 'r', encoding='utf-8') as fasta:
        total_unique_chains = 0  # Use as a flag to test if more than one unique chain
        unique_chains       = []
        for line in fasta:
            if line[0] == '>' and total_unique_chains == 0:
                total_unique_chains += 1
                unique_chains.append(_fasta_helper(line))
                continue
            elif line[0] == '>' and total_unique_chains != 0:
                unique_chains.append(_fasta_helper(line))
    return unique_chains

pdb_diff_chain  = '2p16' # Chain ...
pdb_diff_chains = '1d2z' # Chains ...
pdb_same_chains = '6won' # Chains ...

download_p_diff  = rcsb_downloader(pdb_diff_chain, '.fasta')
download_p_diffs = rcsb_downloader(pdb_diff_chains, '.fasta')
download_p_same  = rcsb_downloader(pdb_same_chains, '.fasta')

unique_s  = fasta_parse(download_p_same)
unique_d  = fasta_parse(download_p_diffs)
unique_ds = fasta_parse(download_p_diff)

print(f'First unique chains in {pdb_diff_chains} (2 unique chains): {unique_ds}\n'
      f'First unique chains in {pdb_diff_chain} (2 unique chains): {unique_d}\n'
      f'First unique chains in {pdb_same_chains} (1 unique chain):  {unique_s}')

```

---
Tags: sequence-alignment, proteins, python, sequence
---


# Question
Title: How many protein sequences do we usually have for each species?
My question is exactly the one in the title. I add a particular example for more clarity. If I want to study the sequence of protein IkBα, do I find only one sequence or is there a database that contains many sequences from many human host species ?

EDIT: The analysis that I am going to do is to use NLP models to apply it to protein sequences to extrapolate, for instance, protein function knowledge.

Here is a reference of a research paper Learning the protein language: Evolution, structure, and function that considers this application by taking into account the problem of low number of sequences data.

So I think I should have at least 10 protein sequences.

Thank you in advance.

# Answer
> 1 votes
Just to add to Maximilian Press's answer. I happen to know about human immunology and I think Maximilian Press has a good idea.

Genetic studies and immunological studies *per se* in cellular immunology tend to be restricted within humans. There are limited exceptions to this namely chimpanzees and humanised mammals, most notably humanised mice. In cellular immunology these mammals retain immunologically distinct behaviour but are sufficiently similar to be useful models.

The protein in question is a B-cell repressor, thats more complicated ...

However, using trans-species homology of the protein of interest for

* chimpanzees
* humanised mice
* regular mice

.. represents a good diversity of proteins to use what appears to be a deep-learning function prediction if this was a traditional cellular immunological study. Any immunologist looking at the biological rationale would recognise the model systems deployed.

**Method: Blast** Generating diversity for given protein is easy, except for humans. Just use online Blast and restrict the search to the species of interest and 20000 hits (?I forget the upper limit) is the maximum and select for more diverse genes. This would work for chimps and (possibly) mice. Thus blast a human gene and restrict the output to the species of interest. Human gene diversity - this much bigger than an online blast could be targeted in accordance with human genetic diversity. The cellular immune system for humans does not conform to background diversity *in general*. Some geographic human populations do have a HLA predominance but there are no absolute association.

The final option is local blasting the human database for that gene, that is doable and objectively used to derive representative diversity - for this purpose is complicated and would advise against it.

# Answer
> 2 votes
I'm going to start by assuming that you don't actually want 10 protein sequences from the same species, even though you write this in your comment on M\__'s (deleted) answer. This will not be a very informative sample of variation in a given protein, except for extreme examples like antibodies.

To get a good sample of variation you almost certainly want to look across species. The paper that you link to does this by using the UniRef database, they then cluster proteins at the 90% similarity threshold. That is well beyond the variation that you'd expect within a species. See their section "Training datasets":

> We train our masked language models on a large corpus of protein sequences, UniRef90 (Suzek et al., 2007), retrieved in July 2018. This dataset contains 76,215,872 protein sequences filtered to 90% sequence identity.

You could use an approach like that in the paper to get similar-ish clusters of proteins whose variation you can then analyze. You can use this, but I am somewhat suspicious of that, without reading the paper in too much depth, because we know rather a lot about orthology/paralogy and how it affects function in these cases. Especially in these very elaborate not very introspectable models I worry about what they will do with inadequately curated data.

So instead I'd probably want to first identify the orthology groups (gene families) of interest. These are generally manually curated groupings of proteins.

Once you have the orthology groups that you care about, you can use a database such as KEGG or OMA to identify the specific orthologs (genes) that belong to the group. Then you will have your dataset.

---
Tags: proteins, sequence-analysis, public-databases, public-dataset, deep-learning
---


# Question
Title: Mapping single cell data with annotation
Entry level Bioinformatician here. I am trying to map my single cell data with standard annotation file or assign my Seurat cluster to cell types. I have used the singleR and celldex package to map it but my problem is the single cell only consist of endothelial cells. I haven't found paper that only consist of rats endothelial cells. The single cell data is for rats. I am stuck. What is the best thing that could be done here?

# Answer
> 3 votes
It seems like you can annotate all your cells as endothelial cells, problem solved ;) But seriously, not clear what your expected cell types are (subtypes of endothelial cells perhaps) or the tissue type. If you can't find a suitable reference dataset for the reference mapping approach, then you can always resort to the 'old fashioned' way of unsupervised clustering (e.g leiden) plus checking marker gene expression to annotate the clusters.

An alternative approach, could be finding a suitable human/mouse reference and then converting the gene ids to rat (ideally via 1:1 ortholog mapping) but this could get messy, so I suggest the former simpler way, at least as a starting point

---
Tags: scrnaseq, clustering
---


# Question
Title: Variant Filtering Using AnthonyOzerov/variant-filtering Github Repository And Python Errors
I am using the following Github repository that I cloned. One apparent issue is that it is using an old/defunct link for the phenotype\_to\_genes.txt that is no longer active and after several tries I found what seems to be the correct link. However, I cannot seem to get it to work as it seems to be throwing up an error that makes no sense to me when I try and run it via a bash script:

**Variable I defined in advance:**

* `pedigree_file` has the path to the pedigree file with 4 columns (Family ID, individual\_ID,Status \[Child,Father,Mother,Sibling\], Sex \[Male or Female\] and Phenotype \[Affected or Unaffected\]
* `Mendelian` a path to a file created via Mendelian filtering
* `output_dir` path to the output directory
* `phenotype_path` a path to a file with 2 columns group/family and then HPO
* `VC_Python_directory` path to the variant filtering directory
* `Map_File` path to the phenotype\_to\_genes.txt
* `family_names` a list of family/group names

**The bash script:**

```
cd "$VC_Python_directory"

IFS=',' read -ra family_array <<< "$family_names"
for family_name in "${family_array[@]}"; do

    # Concatenate the file path
    file_path="$output_dir/$family_name.txt"

    # Echo the family name and file path
    echo "Family Name: ${family_name}"
echo "File Path: ${file_path}"

    # Run the Python script with the correct arguments
    python main.py -p "$pedigree_file" --data "$Mendelian" --output "$file_path" --family "$family_name" -ph "$phenotype_path" --mapfile "$Map_File"
done

```

**The error:**

```
Family Name: Sample_1
File Path: /path/to/output_directory/Sample_1.txt
Path/to/variant-filtering/utils.py:77: DtypeWarning: Columns (2) have mixed types. Specify dtype option on import or set low_memory=False.
  phen_to_genes = pd.read_csv(mapfile, sep = '\t', header = None, comment = '#')
Getting relevant genes for family phenotypes...
Traceback (most recent call last):
  File “path/to/variant-filtering/main.py", line 26, in <module>
    load_phen(families, args.phenfile, args.mapfile)
  File “path/to/variant-filtering/utils.py", line 80, in load_phen
    phen_to_genes.columns = ["HPO-id", "HPO label", "gene-id",
  File “path/to/python3.10/site-packages/pandas/core/generic.py", line 6002, in __setattr__
    return object.__setattr__(self, name, value)
  File "pandas/_libs/properties.pyx", line 69, in pandas._libs.properties.AxisProperty.__set__
  File “/path/to/generic.py", line 730, in _set_axis
    self._mgr.set_axis(axis, labels)
  File “/path/to/managers.py", line 225, in set_axis
    self._validate_set_axis(axis, new_labels)
  File “/path/to/base.py", line 70, in _validate_set_axis
    raise ValueError(
ValueError: Length mismatch: Expected axis has 5 elements, new values have 7 elements

```

I have checked over the pedigree file and map file and phenotype path files but I am really unsure of how to go about this pipeline as I have never seen this kind of error before. I am not sure I can ask this here but does anyone perchance have an inkling as to what is causing the rather peculiar error?

# Answer
> 1 votes
Apparently the format of the `phenotype_to_genes.txt` File changed, while at the time when this pipeline was published the file consisted of 7 columns, it only has 5 now. That's what the error message tells you. It's best to file an issue at the github repository to get in touch with the authors. 1) so that other potential users are also aware of this issue and 2) the authors might be willing to fix this issue/help you get an updated version of the pipeline or an older version of the phenotype-file.

---
Tags: python, bash, github
---


# Question
Title: How to handle technical replicates in cellranger count?
I'm working with cellranger count to annotate some publicly available fastq files (SRX7888073). For each sample, 4 runs are available, which I assume represent technical replicates.

How do you recommend to handle them in cellranger? Are they 4 different samples (i.e. 4 separate runs), or 1 sample with 4 lanes?

Thank you !

# Answer
> 3 votes
It depends what you mean by "technical replicates". If it is the exact same library prep spread out over two lanes, or two flowcells, they must be combined at the fastq stage. Cellranger will naturally do this with fastqs that keep their original Illumina-assigned names, but just have the files split by lane.

If one file has a read from a gene with a particular cell barcode and UMI, and another read in a separate file but **the same library prep** has a read from the same gene, with the same barcode and same UMI, those two must **not** be both counted. So if you combine the counts, that's wrong.

A NextSeq naturally has 4 lanes, so that's almost certainly what you are looking at. The read names would confirm that, unless they've been altered. Cat them all together, or alter their names so that they match the original Illumina specs with lanes 1-4. You might have to do that anyway, cellranger used to be very picky about that, and might still be.

# Answer
> 2 votes
If they are truly technical replicates in the most specific sense (e.g. from the same extraction and library prep, just split up for sequencing), then the counts for each identified cell can be combined (i.e. added together). One way to do this is to concatenate the files prior to running any data processing scripts.

If there is anything different about the runs, then they should be treated as non-technical replicates and properly integrated:

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

We had a situation with Rhapsody single cell data where a library was split into three on three different cell barcode plates, and was showing a small amount of library preparation bias that influenced clustering. Carrying out an integration of the libraries helped to reduce this bias.

# Answer
> 0 votes
After some digging, it seems like each of the 4 runs has I1, R1 and R2 files - they're replicates.

You can either

1. run them separately and then perform `cellranger aggr`; or
2. link to all 12 files using soft-links that follow a proper naming convention (the files need to have the same prefix, lane numbers need to be different per replicate), ensure all softlinks are in the same directory, then use the directory with `cellranger count`.

I'd recommend the former approach as it will probably save a little compute time, is parallelize-able and will be easier to debug/drill down into.

---
Tags: scrnaseq, cellranger
---


title: How do I fix this limma line?

text:
# Question
Title: How do I fix this limma line?
I'm trying to conduct `lmFit` on my `fread::table` matrix:

```
fit <- lmFit(matchedgeneTPM)[,-c(1,2)]

```

It's giving me the error that expression object should be numeric and that there are 2 non-numeric columns (because the first two columns are the labels of the genes and gene ID's I'm working with). I thought the last portion of my code would allow it to skip over the first two columns but apparently it's not the right way to do it. Can someone show me how to format the code correctly, please?

# Answer
> 2 votes
My guess is you need to index `matchedgeneTPM` before you run `lmFit()`, not after:

```
fit <- lmFit(matchedgeneTPM[,-c(1:2)])

```

---
Tags: r, rna, limma
---


# Question
Title: how to search and calculate occurence of certain word (like TAA) in fastafile
```
>003197.2:c1179912-1178219_reverse_complement
AGGCTTTTAAAAGCCTCCTGGGTTTTTAGTGAAGCTGAGTGATAGAAGCTCTGTATTTGCATAGCGTTTT
GTCGGGCGCTTTAGCCGGCGCTTTGCCCTGGCCTGAGAGAATCTGCATTCCGTTGTATAAGGTTTTTTGT
CCTTCTGATGCTGTAGATAATTCCCCCATGTCGCGCCGGGTTCTCGCAGGACAATAATTTCCGGGCGAGC
AGAACGGTTGCCGCGACGGTAAGAAGATCGCGCTGTAAGTTATAGAGGTTATGCAGCGAGTGGTTAGACG
CCGGTCCGCTTTAACTTTGGCTAACTCCATTTGGTTTGCCATTGACGTTAGAACCGGGTCTTGTTTACCC
GTGCTGCAATAAGTTCGATAAGATTTTTCTTCAATGCTTTTGCCGCGGCTTCTTCTTGTTTTGTTGCTGG
GCATCTCTAAACGCTACTGCCGCAAAGCGATGAGCTTCTTTTGCAGGTAAGCCATCCTGCTGCTGAGTGC
GTGATGCCCGTTATGCGTGAGTGTATTTTTTATGGTTTGCCAGGGCTGGTTATTAAGCTGCTTGACCTGA
CATAAGCACTGGGAAAGATATCTTTTGCGCCGATTTTCATCTCTGCTGCAGGGAGCTGCGTGTTGGTATA
ACCGTGGACATCCACAAATTATTGGCGTGATGAATATTCTTGGTATCCCAACTGCATACGCCCTTTCCCT
ATAGGGGGAAAGCACGCCATGACGTATCCCGCAAAAAAGCGTTTTATCTTTACCGTCCTCATGCACACTC
GTGCCGCAGTTAATACTTCTTTGGCTTTGTTTTCAGCGCCGACGTGACGCAGAAGCGGATCTTTTTCATG
ACGGATACCAGTTTCAGGCTTACCGCCTCGCCCGCTAAGGCTTTGTTAAGCAACTCAGGTTTACTAAAAA
CCATGCGCGCATTTGGTCCTCGACCATCGTTCCCTCTTTGCCGAAAATATTCGACGCGGTGAGTAACCCG
CCGTCTGTAGATCGCCATCTTTATTGCGGATTTTTAAATGAATCATTTTTCCCGGCTGGGTCAACGATTG
CCAAAGCCGAGCTTGAGCGCCAGCTCATTAACACCCACATTAAATGCGGCGACGTCCGGTTTTATTTTTA
TTCAGGGCGTAAATCATTGCCTAATAACTGATGTAGCGCCTCGGCATTATAGCTATCCGATGCCTTAAGG
ATGTATTGACGACCTCATAATTATCCGGGTATTGCGCCAGCCATTCGCCAACCCAGCCACCTGGTCTGGC
AGTTTATAGGGTTCGCCGCCATCTTTATGATGTTGGTTATTTTTCCATATATCCTTAATCTGGCGCGCTA
TTTGCCGCTTTTACAATTCCAGGCGGGTACCGCGTCAATTTCATGGGCTAACATGGCAAGGCGTTGTGCG
TTAACATATGGGTCTGATGTAAGGAAATGATCTCTCGCTTGATTTCTGAATCCATCATCCCTGTACGATC
TTACCGCTATTCAGTAATACTTTTTGGAAAATTTTCTGTCCACCGCTATCCGGAAGACTACCAGGCGCAC
CACCTCTGGCGATAAATTTTTCATTACTTTGTTTCCCGCCCCGCCCGTATTTTGTTTCTGAATCTCCAGG
ATGAAGAAATGCCTTTTACTGACTGCCAAATATTTTCATCCCCAACTCGTTTTTGATAGGAAAGATTGAG
TCAAGATGTGATTA
>003197.2:c3026750-3024686_reverse_complement
TTTGATCTCGGTCTGCATACCTGGCATTATGACGGGGGGCTGAGTCCTTACACTTGTAACCATTATTAAT
CCCCTGACAGCGTCGCTGTGGCACTTTCTCTGACTGCGGAAAGATTCGCCGCAAGATTCGTGGCCTGCGT
TTAAACAATGCATCCAGACTCGCCTGTTTGATCAGCGCGGGAAAATCTTCCAGTTGCGGGCCTTTAATTT
CGCTTTCTTACCGGCGACATTATTTGAATTGGTAAAAACTTCTTTTAACGCCTCAGCGTCTTTCCCGCAT
GCGTCTTCGCCTCAGGAGAATCACTGGTGGCGTTGAGCGCTGAACGAAAGAGCCCGGCAAACTCCATTAT
AAAGCTGACTTTTCTTTCAGGCCATCTTTTATGATTTGCGCGGTATACTCTGCCCCCACCTTCATCAGTA
ATCTAACAGCTGCTTTTCCAGGTTTTCCAGCTTCGCTTCCGCTTTCTTTGTTTCTGGCAGCCATGGCCCA
GCATCACCTCCTGCATAACAAGATTCGTCGATAATTTACTTAATTCACCGCCAGTGTTATTTTTGATAAT
GTAAGATTGCTGCGGGTTAACGGATCTAACGTACAGCCAAAGTTATGTTCAATGCAGCTGGCAATATAGG
TTGCGTTAATTTTTGGTGGCACATATCAAGCGCCTCAACAGCCTTCGCCGCCGCTTTGTCAACAAGGTGC
GTAGCAGCAGGGGTATCAACGTTTGCATTTCAAGGTGCCGGGCTTCCCGTCCTACGCTGGTACCCTGCTC
GGTCCAGCCGCAGGCTCCGGAATTTTAGGGTCAGGCAGTTTATCTGCAGGAATTTGTGCAAAAACATTAC
GTTACTGCTATCAATATTAATATTGATATTAATACCCGTAGGTTCTGCTTTTTTCCCACCATCAGGTACT
CGATATGGCGCTGGCTATTATCGACATGGCTTCGGCTATTGTTAATATGCTTACTGTTATCCACGCTATG
CGTTGGCTGTTATCAATAAATGTTCGGCTATTATCAATCGTCTTCCGGCTATTATCATGGTTGCTATTAT
AGTCGAATCCACACGCGAATGACTATGACTTACGTTGCTGCTGTTAGTCGAATGGTGTGACTCGCCATTG
CAATTTTACCCGCGATGCCATGGTCAATTGCCCCTGTGCTGGCGCTGTGTGCCGTCTCCGTTTGATGCGT
ACCACCTTTCCATCTTTACTTTCGGATGAAGCGTTGGTCACAGCCTCTGCCGTCGCATGAGCCGTGACGT
GCTCTTGCTGGTTACGCCATCGACTTCATCAAAAGAGGTTGTTTCACCCGTAGTGCCCTTTTCTGACGTG
CTGCCGTCTGGCTTTGCTGTTTATTATCATCAACGCCGTGCACCGTGGCCTGTACCGGCTTGCCGATAAT
TGTAAGGTATCGGTTTTAACATTCTCGACGCCAGCCAGTTGACTGGCAAGCGGCTTCACATTCACAATCT
CTTATCTGTCTGACTATTATCATTATCGGTCGTACCGGCTTTATTACCTGTAATTACCGTCGTGTCTGAT
TATTCGTAACGCTCGGTACTGTTGAGAGGAAACTATTTTGCTTAAACTTCAAACCAGAAAATGGGCCAGT
GTATCTGGCTCCAGGGCTTTACGTATCACACCCAAAAACGTCTCACGGGCGTCGAAATGCATTGAATGCA
TATCGTGTCATTGGGTAATATTTCCGCACGCAAGCCGTCAAACGCTCTGCGTACAGGGAACGGTGTGGAG
TCTCCATTTTCGCTTTCAATTCTGGGTGCTTGTCAATATCGCTACATTGCGCTTTCAACGCTGCACTCTT
CTGGCAAACTGCAGTGCAATCTCAGCCAGTTTTTCTTTTTGTGGATGATGAGTAATGACCTCTTTGAGAG
GATAAGGCCGTCCAGCACATTACTTAACACATAGTCAGTTTCCCCTTTTAGCCTGGTCAGCCCCGCCTCT
TAACCCGATCCACGCCAGGTTTATTCAGGTATGACTCGTAAGCCGGGCCAGCTCGCCAGCTACCGTCTCC
TTAACGCTGCATGTGCAAGCCATCAACGGTAGTAA
>003197.2:c3010816-3009905_reverse_complement
TATTACGAGTCGTAGGACTTAGCATACTTTTCCCTCCACATGATAGCTCCTGCACCGAAAATATCATCTT
ACAATACAGGTCAATGCTTCGGTGTTTAGCTCCCCGCTGACGTCAGGCTGAGGCGATAAACTCGCCCCCA
TTCGAGATCGGTTTCCTCGTAACTGATATGAATCCAGCTGCCATCTATAATTTCACTTTCCAGACGCTCA
CTGACATAACAAATTTAAGATTTAACTCTGGATACTTCTTATTGGCCTGTGCAACAAGAAAAGGCATCAT
TTTACGAGAAATCGGGCAGATTCAACGCCTTCCATTCTGGTCTCCTTTATAGAGGAAACAAGCTCATGGA
AATCACCGATGTCTTTCCGTCCATAACTCTAAAATCGACGACTGAAATATGGATACCTGAACTTCCCATA
CAAGAGCTGCTTTGGTCCTCAACGCCAGTAAAGCAGGTCCAAAAGCGCTACACGCTGCTGGTTCGAACAA
CATTCAGAAGAGCTTCGTTGAATGTCCAGCTCGACCATAGCAAAATAACAATCAGGCAGTTGTTCACGTT
TTCATGAATTTTTACTAAGTTCATAAATTCAAGCTGAAGTTTTTTGGCGAGCGCCAGGCTAAAAATACCG
GCAGATAGCGGTCTGCTTTATCGGACGGGAGAAAAGGTTCTTCACCACACAGACGTTCACAAATATTATC
GCCGCCGGATTGGCCTCCACATATTCATTTAATCGTTGTGCGCCTTGAGTATGTTTGTAAAAACTCACCG
TAGCATAACGGCATTGTTATCGAATCGCTCATAAAGCGTTTCATTCTTTTTGTTTACTATACTGCTTCCC
ACGGTTTAAGTAAAGACTTATATTCAGCTATCCTTTTTTTATGAGCGGATATAGAGAGTTTTTTATCGTT
TT
>003277.2:c27903-27178_reverse_complement
CCCCATCATAAGCAGCTACAATATTCAAAGGAGGGATGTTTAGATTTAGATTAGGCCTATTTATGGGCAT
ATCTTCCGCATTTGGTTGTGCAAGGAGTTGAAATTATTTTTCAGGTGCTTATTTGTAGATGGTATTTCCG
AAAGAAACCACTCTGGCGCATCCCTGAATAGTCAGGCACATCCAGCGCCTCTTTAAAGTGGGATACCGGC
GTGATTCTCGACGAGCATGGATGAAAACATCGTAACCATGGTTATTCAGCTGAAAACCTTGGCTCATAGC
GCTTGTGGCACCATATCCCTGAGCACACTGATGTGGAACTTGTCACCGGCAAATTTGCCCTGAGACTGAG
CATATCGGTCACTTTCCACTTATCTACCGGACTGTCCTCTGAAAACAGCAATCCGGACAGCGCTTGAAAT
CCTGGTCTGGTTTTATATACAACGTGAACTGAGCGCCCAGGCTAACACGGGCTTGTTGAACGACCTTCTC
AGTCTGGATTCCAGACACTCTATAAATTGCCGTGTCTTGTGGAGAAACGACGCACTGTACTGCGAATTTT
GAGATATTTCCAATTTTCAGGATGAACGTCTGACTCAGGACACTGTCCTGAAATAACCCCATTTTCGGAT
GTTCCTCACGTAAAGCCTGTCTCTGCATTTCGCCACCATCACGCCCACTTCGTAGTTCATTACGATAACT
TTACTCTGTCATCAAACGATAAAACG

```

# Answer
> 1 votes
Since you tag python, I suggest a simple Biopython solution:

```
from Bio import SeqIO

print(sum(s.seq.count("TAA") for s in SeqIO.parse("file.fa", "fasta")))
# 97 on your input file

```

If you wanted an overlapping sequence count, Biopython has a method for this too:

```
print(sum(s.seq.count_overlap("TAA") for s in SeqIO.parse("file.fa", "fasta")))

```

Here they are the same, but for a string like "TAT" they could be different:

```
>>> from Bio.Seq import Seq
>>> Seq("TATAT").count("TAT")
1
>>> Seq("TATAT").count_overlap("TAT")
2

```

---
Tags: python, fasta, multi-fasta
---


title: How do I calculate DEG cutoffs?

text:
# Question
Title: How do I calculate DEG cutoffs?
So for context, I have a set of TPM values for multiple genes for different samples, and I need to calculate the differential expression for RNA-seq values. I've been using this website as a guide: https://raw.githubusercontent.com/ucdavis-bioinformatics-training/2019\_March\_UCSF\_mRNAseq\_Workshop/master/differential\_expression/DE\_Analysis.Rmd. The overall aim of this calculation is so that I can use PANDA to map TM motifs/combine gene expression data to create networks representing interactions between transcription factors and genes, and GSEA to analyze genes ranked by fold change or differential expression p-value.

Note that I am using log2(TPM+1) values. How do I calculate what the filtering cutoff should be? Note that I already excluded all genes that had any values of 0 for any of the samples, and I will use the Benjamini-Hochberg method to adjust p-values for multiple testing afterwards.

# Answer
tl:dr use limma (not voom, that's for raw counts) and refer to previous questions:

The problem you have is not so much the threshold but the method to infer any metric to be used for this threshold. The arguably most common methods in R for DEG analysis are DESeq2, edgeR and limma, based on number of citations. The former two require raw counts which you apparently do not have, but if you do then use this with any of the mentioned tools.

limma is what is left and could be used here with TPM. Note that TPM has the problem that counts are divided by a length factor. A brief background: Any statistical test aims to decide whether the given variance between two groups (basically the logFC) is greater than expected by chance, which is the variance within group. In RNA-seq the variance is not equal across the range of counts (low, middle, high, for simplicity). On log2 count scale variance is greater for genes with low(er) counts than for higher counts. Usually, using raw counts the method (DESeq2, edgeR, limma, and others) learn from the data what is the expected variance at the given average expression level, and incorporate that into the model to make adequate inference even at small sample size. That is by the way where the common assumption "most genes do not change" comes from. Assuming that most genes are not DE, then the mean-variance trend represents the expected variance at the given average count range, and this can then be used to help decide whether a gene is DE or not. With TPM you basically throw this out of the window because the actual counts (as they come from the sequencer) are manipulated by dividing by length and the entire mean-variance model gets distorted. Another factor is that TPM, which corrects by total library size can fail to properly normalize counts, I wrote a summary here:

**Long story short: TPM is not good for DEG.**

If it's all you have then use what is called "limma-trend", see the limma manual at Bioconductor, which will directly model the log2-TPMs and give some logFC and FDR. You can then use standard cutoffs, such as FDR \< 0.05 & logFC \>/\< 1.5. These cutoffs are not perfect, but common, and given that TPMs are bad anyway you can and should stick to standard cutoffs as much as possible, rather than infering your own ones. This at least gives some standardization to this suboptimal analysis.
> 1 votes

# Answer
Two solutions - the underlying issue is called Type 1 error (below).

1. If you are performing studies on humans then the Benjamini-Hochberg method is without question the method of choice, because this includes the *false detection* rate for EdgeR and DeSeq2. Thats been independently calculated.

Thus, the paper that caused a stir was Li *et al* (2022) who wrote,

> Particularly, the actual FDRs of DESeq2 and edgeR sometimes exceed 20% when the target FDR is 5%.

***Benjamini-Hochberg (BH) method solves that for humans (EdgeR and DeSeq2)***, because the problem with the BH method is you need a false discovery rate, Li *et al* 2022, calculated that its about 20% ... so for EdgeR and DeSeq2 the problem is solved. The BH tests will simply use the estimate for the FDR by Li *et al* (2022) instead of the FDR by EdgeR and DeSeq2. Cool? It is :-) it wasn't what Li *et al* (2022) intended but it is the solution to a difficult conundrum (concerning a statistical issue or parametrics vs. non-parametrics). **So yeah it's cool** (needs count data).

To perform BH you simply rank the observed probabilities (e.g. from Limma), and there are significant if each is below (or equal to) the BH rank. The BH rank is the total number of tests as denominator, the rank position as numerator and this is multiplied by the false detection rate as a frequency.

*HOWEVER*, you are using Limma and I don't how the model works. I know the underlying stats behind EdgeR, DeSeq2 and NOISeq.

---

Okay so the false discovery rate is not known ... in the way Li *et al* (2022) calculated it. What now? Limma is likely to calculate an FDR, you could use that via BH. My system is so far from humans there's no way Li *et al* (2022) FDR is relevant.

2. Bonferoni's approach - I don't know how mainstream this is, but its what I use (I've a background in classical stats), because I've no idea what the true FDR is. This approach can be exemplified via Volcano plots is to visualise the Type 1 error issues.

Here's a volcano plot which permits online viewing via links obtained here. I can't change the size BTW (Bioconductor and loads more BTW don't permit this),

This is using a flat probability - which in the picture is very low - it's the horizontal threshold (*y*-axis). So the P value is one parameter approach to significance ....

The other parameter is the log2 threshold ...

Here's a two-parameter approach obtained from GeneKitR's site, who do permit images to be link viewed (thank you genekitr).

Here the authors have used log2 as the threshold (*x*-axis). The threshold is probably log2 \> 1 and log2 \< -1, which is reasonable. They have set the P-value threshold to be really low, looks to be about 5%. The plot reveals the effective P-value is 1E-2.5, which is low, but nothing really starts happening until 1E-3, which ... okay. Although this is a two-parameter approach the guys really just relied on log2.

**Type 1 error and a 1-parameter method threshold**

The massive problem is the shear number of tests performed and this will result in significant values PURELY BY CHANCE ...artefacts of critical probability. BH is one approach, Bonferoni's is an alternative way in.

Put simply for your data to calculate the Bonferoni threshold if a 5% critical probability is sought then its ...

> 0.05 / 19000 = 2.6e-6

Which is okay. So any single result \<2.6E-6 is significant. This okay, but you will lose results so its vulnerable to Type 2 error. Why do you lose stuff ...

**Type 2 error and 1-parameter method threshold**

The biological question is e.g. is 0.1log2 biological sensible? What does it actually mean? Honestly I don't have a clue ... and the problem with a 1-parameter approach, it's very conservative but infers weird stuff. An expression difference of 10% isn't really serious, but thats what it's measuring. ... because its measuring this it means the probability via Bonferoni's is set high - *for meaningless stuff* and you can be missing important shifts in log2, e.g. \>1 because they are falling under the high threshold.

**Type 1 error and a 2-parameter method thresholds**

This is where it all starts getting interesting. As soon as a +/- log2 threshold is set the number of total tests will drop A LOT. That means the P-value threshold raises. For example, if selecting \> (or \<) +/- 1log2 as used in the above example, supposing the number of observations above/below the +/- 1log2 was 5000. The P-value now raises via Bonferoni, e.g.

> 0.05/5000 = 1e-05

That looks better there'll be more results. Supposing we set a threshold of +/- 1.5log2 and obtained 500 results ... then the new Bonferoni derived threshold is,

> 0.05/500 = 1e-04 (0.0001)

That looks fabulous giving loads of significant results ... Wait isn't this all getting a bit crazy, because it depends on the biological assumption of what a reasonable log2 is. Yep, indeed, agreed.

**Type 1 error and a 2-parameter method thresholds: weaknesses**

The issue is that a reasonable log2 is sensible, because there's a point where log2 expression differences becoming utterly meaningless. However there's a point up to -/+ 2.5 log2 (which is very high) where you could choose anything - that shifts the critical probability and that shifts how many results are positive.

So how do you ...

**Set a log2 within a 2-parameter threshold?**

I've just solved that. I was really happy. The practical answer is it's complicated. My contractual basis is I can't reveal that because public disclosure isn't permitted without top-level approval.

The approach is to use an independent method to calculation DEG (if you check my posts you'll see I do AI). You can then perform parameter tuning of DEG against the independent method and statistically assess the association. In my case I was surprised by the answer - intuitively I wouldn't have gone there, but my system ain't humans.

**Conclusion**

In your case,

* simply check the literature to see what other log2 thresholds are in use, use that and adjust the P-value threshold accordingly as described above and cite that work.
* look at the biological reality of the results at different log2 thresholds. For example, an operon should behave like an operon. If the log2 threshold is smashing-up that operon - that doesn't make sense and adjust log2 threshold by lowering.
* alternatively use BH, this however assumes Limma's FDR calculation is correct and the problem is there no independent study examining that.

Either of the above 3 approaches will be fine your context. This is not cutting-edge statistical methodology it simply needs to be reasonable.

**Ultra-conservative parameters and Type 2 error**

To avoid all the complexities some investigators simply use \> (or \<) +/- 2.5log2 and P-value = 1e-6. These values in ANY system are extreme when used together, but there is no chance of Type 1 error. The danger is Type 2 error , where significant values are either omitted or considered non-significant. This could happen here from two fronts:

* If this results e.g. 100 observations, and DEGs with P-values \< 5e-4 but \> 1e-6 are omitted
* DEGs \< (or \>) +/- 2.5 log2 and above the lower threshold are omitted.

However, it will depend on the question.

**Note** Much of the theory in *point 2*, would also apply to the BH-method. However, I've worked this theory exclusively for Bonferoni's correction and I've not given much thought to BH. This is because it's not relevant to my purpose.

---

*See comments below your post about your TPM approach*. This post is cut-off criteria.

---

Just to finish this post, here's a volcano with thanks to Britannica for permitting the link-view. I think the idea of a volcano plot refers to the "plume" of ash, dust and occasionally lava (which some of us have real-world experience).
> 0 votes

---
Tags: differential-expression
---


title: ZIP the fasterq-dump files in loop

text:
# Question
Title: ZIP the fasterq-dump files in loop
I want to zip the files generated by fasterq-dump simultaneously as they are generated. The code that I tried is this.

```
cat SRR_Acc_List.txt | while read i
do  fasterq-dump ${i} --split-spot -v --skip-technical  -e 15| gzip -k  ${i}> ${i}.gz
done

```

# Answer
Haven't tried, but I believe you may add `-Z` to fasterq-dump and remove `-k ${i}` from gzip. In addition, I would recommend the following:

```
cat SRR_Acc_List.txt | xargs -i echo fasterq-dump {} --split-spot -Z -v --skip-technical -e 15 \| gzip > {}.gz | sh

```

It is simpler and more flexible. For example, replacing `sh` with `parallel` would allow you to download in parallel.
> 1 votes

# Answer
I tried this and it worked

```
cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-spot -v --skip-technical  -e 20 ; for j in i ; do   gzip ${i}.fastq  ; done; done

```
> 1 votes

---
Tags: ncbi, sra
---


# Question
Title: How do I correctly format my limma/eBayes code?
Here is part of my data table, which I will call "df":

```
    Subject.ID    Sex           V1           V2           V3           V4           V5
1   GTEX-1117F female   4.24944534   0.18307358   1.99276843   0.21785110   0.74777407
2   GTEX-1128S male    3.286585483  0.165944088  1.843983844  0.444455046  0.307311939
3   GTEX-11EMC female   3.31947343   0.27846061   1.46519089   0.33753984   0.91708799
4   GTEX-11GSP female  3.232353768  0.215492230  1.778208576  0.166458372  0.644225472
5   GTEX-11TTK male    2.934705200  0.256406854  1.712815854  0.256165277  0.948301812
6   GTEX-11ZUS male     3.77188558   0.18434378   2.43189035   0.37373172   1.14339342
7   GTEX-12WSD female   4.22032995   0.21933892   0.93900085   0.08687886   1.06901468
8   GTEX-12ZZX female   3.43616184   0.39473363   2.73205207   0.33525457   1.51399598
9   GTEX-1313W male     3.66334462   0.39418485   2.47248777   0.53545580   0.64634702
10  GTEX-131XW female  2.956614288  0.276317993  1.977096712  0.167743280  1.528071165

```

How do I correctly format the following code to account for the kind of dataframe I'm working with? I'm using sex as the factors to be interacted. Here is what I have so far:

```
design <- model.matrix(~ Sex)
fit <- lmFit(df, design)
fit <- eBayes(fit)
topTable(fit)

```

The first line is giving me the error "Error in eval(predvars, data, env) : object 'Sex' not found" (and idk why I'm getting this error since there is a "Sex" column in df). The second line gives me the error "Expression object should be numeric, instead it is a data.frame with n non-numeric columns" (and I know why I'm getting this error; I'm not sure how to handle the error, though).

# Answer
You need to look at how to create a model matrix : https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/model.matrix

In your code, you did not specify the data from which the Sex columns should be taken, or you did not created the Sex variable

Either create `Sex=df$Sex` or specify it in the model.matrix function `model.matrix(~ Sex, data=df)`
> 2 votes

---
Tags: r, limma
---


# Question
Title: Removing missing categories from geom_bar
I have a simple dataframe that includes points where the x values are 0, 1, 2, 3, 4, and 7. When I go to plot the data this is what ggplot2 gives me:

I only want to include the points mentioned above. How do I take out the instance of 5 and 6 in the bar plot?

(I have tried googling the answer and nothing is coming up with my actual need so I figured I'd ask the community which might give an answer sooner - please be kind I know this is bound to have a simple solution)

# Answer
> 4 votes
I was able to get it to work by converting the column in the dataframe to a character using this command:

`df$day <- as.character(df$day)`

and it worked like a charm! Thanks for your advice @haci

> You can give a try by converting the column of numbers to a factor variable.

---
Tags: r, ggplot2
---


title: How to parallelize BEAST

text:
# Question
Title: How to parallelize BEAST
I'm running **BEAST v1.10.4** with **BEAGLE v4.0.0**. My dataset has one large partition and I have access to a machine with many CPUs so I would like to optimize performance by parallelizing the process. Following the documentation and various (quite old and probably outdated) help Q&A pages I tried

```
beast -beagle -beagle_CPU -beagle_instances X -beagle_thread_count Y -threads Z mydata.xml

```

with different values for XYZ but BEAST is ever only using a single CPU regardless.

**How do I make BEAST use on more than one CPU?**

We have a 'server' in-house, that is essentially a big computer with 144 CPUs that is shared by all users. That's why it would be important to limit the max amount of threads with the -threads option. No queuing system, no GPUs. I want to do some tip date sampling which afaik is not available in Beast2.

# Answer
The machine is important because of the introduction of M1 type chips which will run GPU.

From historical use of beast1 the following will work ...

```
beast -overwrite -beagle_auto -seed 2001 mydata.xml > myfile.runtime.log

```

This will parallelise happily within the number of processors available. I would recommend it. This is because it may not parallelise across 144 CPUs. It depends on the machine architecture. Also from memory (mine not RAM) it dynamically readjusts, so even if it was 144 CPUs parallelised when another job lands on a CPU that drops to 143 CPUs for beast1.

Anyway, I use the Cloud so I no longer keep up with computer architecture.

If this floods 144 processors and there's no dynamic readjustment I would suggest,

```
 beast -beagle -beagle_SSE -beagle_CPU -threads 10 -overwrite mydata.xml > myfile.runtime.log

```

I would limit the number of threads to get it working before considering e.g. 40.

However, shifting to GPU is strongly recommended, this will solve any beast2 issues and its very fast. However, I don't know about beast1 for GPU, I've never used it for that.

---

If none of the above work and you've not get any sort of M1 or higher chip (to move to GPU) then you need to look at the error log, especially if you see this error.

> Error 122 parsing the xml input file Cannot create class: ThreadedTreeLikelihood.

---

**Note** although you've got a nice CPU heavy machine, simple fairly modern Mac will do the same operation probably quicker via GPU, at least for beast2, beast1 I dunno

---

Just to note I'm acting in solidarity with an SE strike until an internal matter is resolved about ChatGPT and requested a meeting with staff to explain why. I've formally notified SE. I'll be answering very select questions and doing only essential admin (breaches of CoC and spam). Hopefully all is back to normal shortly.
> 1 votes

---
Tags: phylogenetics, phylogeny, beast
---


# Question
Title: How to identify index sequences for cutadapt for atac-seq
I am just starting on my atac-seq analysis (not very experienced, so my apologies for a stupid question), and I am wanting to trim the adapters of off my paired-end reads.

I am following the galaxy atac-seq protocol that says to enter in the custom 3' adapter sequence. I took the primer sequences from Buenrostro.

```
>Ad1_noMX
AATGATACGGCGACCACCGAGATCTACACTCGTCGGCAGCGTCAGATGTG
>Ad2.7_CTCTCTAC
CAAGCAGAAGACGGCATACGAGATGTAGAGAGGTCTCGTGGGCTCGGAGATGT

```

Following the ATAC protocol, I see that the primer sequences they used are the reverse complement sequences of the end of the provided sequences (ACATCTCCGAGCCCACGAGAC/CACATCTGACGCTGCCGACGA). I took this part and added the reverse complement index primer of 2.7 (GTAGAGAG) to the front.

So for my R1 my sequence is: GTAGAGAGACATCTCCGAGCCCACGAGAC And for my R2 my sequence is: GTAGAGAGCACATCTGACGCTGCCGACGA

However, when I run the command I only get 1.4% read1/2 with adapter. Could anyone suggest where I am going wrong?

Thank you

# Answer
> 2 votes
The sequence to trim for ATAC-seq is `CTGTCTCTTATACACATCT`. It's the Nextera adapter sequence. The sequence is the reverse complement of the 3' ends of your primers, both forward and reverse.

You trim adapters which come before any index sequence so trimming adapter (and by this everything beyond it) is sufficient.

---
Tags: atac-seq, cutadapt, galaxy, indexing
---


# Question
Title: How to download the whole BLAST nt database into a specific folder?
I have successfully downloaded the whole nt BLAST database into the current folder using:

```
wget -b "ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.??.tar.gz"

```

However, I would like to download it in a specific folder, e.g. output/

I have tried:

```
wget -b "ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.??.tar.gz"  -O output/

```

but it does not work. Nothing happens. My 'wget-log' is empty.

I could move to that folder and download it there but I would like to remain in the same working directory.

C.

# Answer
`-O` is to specify the name of the output *file.* You want the `-P` option.

So, from your example

```
wget -b "ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.??.tar.gz" -P output/

```
> 9 votes

# Answer
I recommend using NCBI's tool for this called `update_blastdb.pl`, shipped with the blast suite (I use the conda environment from bioconda).

`cd` to the directory you want nt DB in and run:

```
update_blastdb.pl --decompress nt

```

Pros: if a download fails for whatever reason it is retried, downloads are md5sum checked and re-running in a folder with an existing DB will (as the name suggests) only update the relevant DB portions. More information here!
> 0 votes

---
Tags: blast, data-download
---


# Question
Title: How are the "effect_allele" and "other_allele" bases chosen in a PGS file from the PGS Catalog?
How are the "effect\_allele" and "other\_allele" bases chosen in a PGS file from the PGS Catalog? For example, for PGS002723 (https://www.pgscatalog.org/score/PGS002723/), if I set the "effect\_allele" to "ref", then "ref" is equal to the reference base of hg38 in 98.60% of the cases on chromosome 22. On chromosome 22, there are 205 PGS variants where the "effect\_allele" does not match the reference base of hg38. In all of these 205 cases, the "other\_allele" matches the hg38 base.

Why does PGS Catalog encode the hg38 reference base in 14,477 of cases for chromosome 22 as "effect\_allele" and in 205 of cases as "other\_allele"?

Finally, how can I deal with this encoding if I want to compute the polygenic score on a certain cohort? Can I simply invert the sign of the "effect\_weight" column for these 205 cases where the "other\_allele" matches the reference?

I paste my reproducible code below:

```
suppressMessages(library(BSgenome.Hsapiens.UCSC.hg38))
suppressMessages(library(data.table))
suppressMessages(library(R.utils))
suppressMessages(library(stringr))

setwd('/cluster/home/criccio/PRS/')

source('scripts/get_hg38_base.R')

if (!dir.exists('output/table'))  dir.create('output/table')

PGS_ID <- 'PGS002723'

# Example to ask questions online
system(paste0('wget -P output/ --recursive --no-parent ftp://ftp.ebi.ac.uk/pub/databases/spot/pgs/scores/', PGS_ID, '/ScoringFiles/'))

scoring_file <- data.table::fread(paste0('output/ftp.ebi.ac.uk/pub/databases/spot/pgs/scores/', PGS_ID, '/ScoringFiles/Harmonized/', PGS_ID, '_hmPOS_GRCh38.txt.gz'))
data.table::setnames(scoring_file, c('chr_name', 'hm_pos', 'effect_allele', 'other_allele'), c('chr', 'pos', 'ref', 'alt'))

pgs_variants_per_chrom <- scoring_file[, .N, chr]
data.table::setnames(pgs_variants_per_chrom, 'N', 'nb_pgs_in_scoring_file')

chromosome_number <- 22

scoring_file[chr == chromosome_number, ref_hg38 := sapply(scoring_file[chr == chromosome_number, pos], get_hg38_chr22_base)]
scoring_file[chr == chromosome_number, ref_hg38 := toupper(ref_hg38)]
scoring_file[chr == chromosome_number, table(ref == ref_hg38)] # 14477 TRUE and 205 FALSE
scoring_file[chr == chromosome_number, 100 * mean(ref == ref_hg38, na.rm = TRUE)] # 98.60%

scoring_file_mismatch <- scoring_file[chr == chromosome_number & ref != ref_hg38]
scoring_file_mismatch[chr == chromosome_number, mean(alt == ref_hg38)] # 100% match
data.table::fwrite(scoring_file_mismatch, paste0('output/table/', PGS_ID, '_chr_pos_ref_mismatching_hg38_chr', chromosome_number, '.csv'), sep = ',')

```

get\_hg38\_base.R

```
library(BSgenome.Hsapiens.UCSC.hg38)
Hsapiens_chr22 <- strsplit(as.character(Hsapiens$chr22), '')[[1]]
get_hg38_chr22_base <- function(pos) Hsapiens_chr22[pos]

```

# Answer
Effect allele is almost always the reference allele. It theoretically does not have to be, so it's good to not rely on that fact. However, in practice it is typically a safe bet. At first glance, the answer to your question could be this:

The polygenic score from the PGS catalog is using GRCh37. Your genome build is hg38. I expect that some alleles in the score swapped ref/alt positions between builds, resulting in your 205 cases where the effect allele matches the alt allele.

To confirm, could you list a few examples of cases that matched as expected and a few examples of cases that didn't match?

You don't need to deal with this encoding if you want to compute a polygenic score on a cohort. You can match based on the allele itself. (Also, there are existing tools which will calculate polygenic scores for you.) Here is an example (just a snippet, you'd have to iterate through your files, etc. for this to work):

```
            if PGS_row scoring file identifier like rsID or position == vcf_row identifier: # When the SNP is in both files
            effect_allele = PGS_row[2] # Wherever the affect allele is in your scoring file
            effect_weight = float(PGS_row[4])
            dosage (of ref allele) = (depends on your VCF formatting)

            if effect_allele == vcf_row[3]: # Effect allele matches reference allele
                running_total += (2 - dosage) * effect_weight # We assume dosage of alt is 2 minus dosage of ref. So if you have 2 alt alleles, 2-2 = 0 and there will be nothing added for this SNP.
                num_additions += 1
            elif effect_allele == vcf_row[4]: # Effect allele matches alternate allele
                running_total += dosage * effect_weight
                num_additions += 1
            else:
                nomatch = nomatch + 1

# Calculate final output score
final_output_score = running_total / num_additions
percent_matched = (num_additions)/((nomatch)+(num_additions))

```

You shouldn't invert the sign of the "effect\_weight" column for these 205 cases where the "other\_allele" matches the reference. Not having the effect allele doesn't add the opposite beta to the score, it simply doesn't change it.

Try this:

1. Do you still see those 205 cases if you use a genome build matching the PGS file? (I.e., GRCh37.)
2. Do you still see those 205 cases if you use the LiftOver PGS scoring file with hg38 coordinates? (Download it here.)
> 2 votes

---
Tags: genomics, gwas
---


title: How do I lift GWAS results to hg38?

text:
# Question
Title: How do I lift GWAS results to hg38?
There are some questions about lifting between reference builds, e.g. this one. But there doesn't appear to be a question about lifting a GWAS results file to a new reference build (except off-site). GWAS results files are not VCFs or plink files, so they aren't covered by the usual lifting tools (such as CrossMap or LiftOver). How do I lift a GWAS results file to hg38?

For instance, the Lee et al 2018 GWAS on education has this file:

```
user@computer:/home/user/Lee_2018$ head GWAS_EA.to10K.txt
MarkerName  CHR POS A1  A2  EAF Beta    SE  Pval
rs9859556   3   49455986    G   T   0.6905  -0.02901    0.00151 4.61e-82
rs7623659   3   49414791    T   C   0.3095  0.02899 0.00151 6.05e-82
rs11917431  3   49644012    C   T   0.6973  -0.0292 0.00152 7.40e-82
rs1873625   3   49666964    C   A   0.6973  -0.02916    0.00152 8.10e-82
rs11921590  3   49644193    T   C   0.6973  -0.02916    0.00152 1.17e-81
rs2352974   3   49890613    T   C   0.5 -0.02569    0.00141 6.35e-74
rs7029718   9   23358495    G   A   0.5646  -0.0244 0.00143 7.82e-65
rs7868984   9   23357826    C   T   0.4388  0.02437 0.00143 8.21e-65
rs11793831  9   23362311    T   G   0.4354  0.02418 0.00143 9.34e-64

```

# Answer
> 3 votes
Yes. LiftOver tools, if they work on VCF data (some only take BED as input) work on GWAS summary statistics. You could use https://pypi.org/project/liftover/, for example.

# Answer
> 2 votes
It turned out to be easier than expected. You figure out which dbSNP version they used. dbSNP is the classification system for snps based on the rsid format. Pretty much every GWAS publishes their results with the rsids for the variants, as done above as well (MarkerName). So you download one of the dbSNP files from the official site. In this case, I am using version 142 since this was used by the 1000 genomes dataset.

After this, it is a matter of joining the tables. This can be done easily in R (you need about 15GB free memory), we load the dbSNP database into memory. Since we don't need all the columns, we skip some of them to save memory:

```
#library
library(tidyverse)

#cols
vcf_cols = cols_only(
  `#CHROM` = col_character(),
  `POS` = col_number(),
  `ID` = col_character(),
  `REF` = col_character(),
  `ALT` = col_character()
  )

#this is from the 1kg
#ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/other_mapping_resources/
dbsnp142 = readr::read_tsv(
  "ALL_20141222.dbSNP142_human_GRCh38.snps.vcf"),
  comment = "##",
  col_types = vcf_cols
  )

GWAS_EA = read_tsv(GWAS_EA.to10K.txt,
                   col_types = cols(chrpos37 = col_character(),
                                    chrpos38 = col_character()))

#make chrpos
dbsnp142 %<>% mutate(
  chrpos38 = str_c((CHROM %>% str_replace("chr", "")), ":", POS)
)

#count overlap
(GWAS_EA$MarkerName %in% dbsnp142$ID) %>% table2()

#mutate and merge
GWAS_EA = GWAS_EA %>%
  mutate(
    chrpos37 = str_c(CHR, ":", POS)
  ) %>%
  left_join(dbsnp142 %>% select(ID, chrpos38), by = c("MarkerName" = "ID"))

#update file on disk
GWAS_EA %>% write_tsv("GWAS_EA.to10K.txt")

```

The two new columns, chrpos37 and chrpos38, have the chr:pos variant names in both versions.

# Answer
> 2 votes
I would advise against using dbSNP rsIDs to perform a liftover as they get updated over time. I recently wrote a software to liftover summary statistics and the documentation includes instructions that show you how to easily generate educational attainment summary statistics for GRCh38

---
Tags: gwas, liftover
---


# Question
Title: Looking for datasets of RNA-Seq (or scRNA-Seq) for L1 larvae C. elegans
I was wondering if anyone know of a repository for C. elegans data sets, specifically in L1 larvae stage.

Are there any papers or GEO-sets for this organism at this stage?

thanks

Assa

# Answer
See WormBase, which has a number of different resources for different nematodes:

https://wormbase.org
> 0 votes

---
Tags: rna-seq, scrnaseq, wormbase
---


# Question
Title: How is the odds ratio of disease risk conferred by a 1-standard deviation increase in PRS calculated?
One standard deviation from the mean is commonly used to calculate a polygenic risk score for GWAS, e.g. human genetic disease.

Why is this a common metric, for example why not 2-SD or 1.96 SD as in the normal distribution and what is the formula?

Here's a figure taken from this paper (10.1186/s13073-020-00742-5) to help guide the explanation:

# Answer
> 1 votes
The "OR/SD" metric is somewhat arbitrary. It's just convenient that it is 1 SD; there is no exact reason. If you wanted to make it OR per 2 SDs, you could just multiply the OR by 2.

Imagine you have a distribution from 1st to 100th percentile of PRS, and you have a linear association between percentile of PRS and probability of disease. How would you quantify this? OR/SD is a quick and easily interpretable way to say how sensitive the PRS is. In other words, how "fast" the OR increases as you increase PRS.

The method for obtaining this metric usually involves converting your raw scores into Z-scores. The slope of your regression where Y is OR and X is Z-score I think would correspond to OR/SD.

---
Tags: gwas
---


# Question
Title: Polygenic Risk score calculation using vcf files
I have about 750 vcf files and I am looking to generate a Polygenic Risk Score for each individual. I am planning to use either PLINK or LDPhred to generate a Risk Score. I have phenotype information of these individuals.

Can anyone verify if the polygenic risk score is calculated using individual vcf files or should I look into combining the vcf files?

# Answer
For the tools you mentioned, you should convert to BED format. Other tools will accept VCF. It should work either way, combined or not, though you may have to adjust your pipeline depending on which you prefer. I'd imagine a combined file would be easier to work with.
> 0 votes

---
Tags: vcf, plink, polygenic-risk-score
---


# Question
Title: Loss of predictive power of polygenic risk score when dataset contains missing variants
I am trying to calculate polygenic risk scores (PRS) scores for a new dataset. This dataset does not have all the variants that the PRS score needs. The PRS score I am interested in has 40 variants, but in my dataset only 20 of those variants are present. What is the correct way of assessing the decreased predictive power of the PRS score due to the missing variants?

I assume it is not as easy as summing the effective weights of those variants I have and comparing to the full sum (all variants needed for the PRS).

Thanks a lot and sorry for newbie question!

# Answer
In a sample with all 40, determine the accuracy of the PRS. Then in the same sample, remove the 20 you are missing in the other set, and then determine the accuracy with the missing variants. This will give you an idea of the level of attenuation. However, you may not have access to such a dataset. In that case, I would recommend imputing your data before calculating a polygenic score. (It's likely that the GWAS the PRS is from used imputation itself.)
> 1 votes

# Answer
Which software package around you using, PRSice? Generally missing values comprising half the total observations are extremely difficult to calculate using regression analysis (or many statistical analyses for that matter.

One way forward is to include (an) example(s) of all 40 variants into your data set and perform the analysis with and without the added data, then compare this against a data set where all 40 variants are observed and your samples are excluded. You might find the missing values have low probability of PRS and this would then explain your results.

Missing analysis in regression is a specialist area of statistics and is not for the faint hearted, so a work-around is the best way forward.
> 0 votes

---
Tags: snp, gwas, polygenic-risk-score
---


# Question
Title: Nextflow - Process has already been used
I am trying to write a nextflow script below:

```
/*
 * pipeline input parameters
 */
params.reads = "/path/to/fq/files/"
params.guide_library = "/path/to/guidelibrary/"
params.outdir = "results"

log.info """\
    DUAL GUIDE CRISPR QC - N F   P I P E L I N E
    ===================================
    guide library: ${params.guide_library}
reads        : ${params.reads}
    outdir       : ${params.outdir}
    """
    .stripIndent(true)

/*
 * define the `index` process that creates a binary index
 * given the transcriptome file
 */

process READCOUNT {
    input:
    path guide_library
    tuple val(sample_name) path(reads)

    output:
    path "$sample_name"

    script:
    """
    exactcounts.py -inputlibrary $guide_library -read1 ${reads[0]} -read2 ${reads[1]} -prefix $sample_name
    """
}

workflow {
    Channel
        .fromFilePairs(params.reads, checkIfExists: true)
        .set { read_pairs_ch }
   count_ch = READCOUNT(params.guide_library,read_pairs_ch)

}

```

But am getting the following error:

```
Process 'READCOUNT' has been already used -- If you need to reuse the same component, include it with a different name or include it in a different workflow context

```

# Answer
> 5 votes
I believe adding a comma after `val(sample_name)` on line `tuple val(sample_name) path(reads)` will solve the issue.

EDIT: Also documented here: https://github.com/nextflow-io/nextflow/discussions/2714

---
Tags: nextflow
---


# Question
Title: Why is my centred log ratio transformation on a distance matrix in R not symmetrical?
I am trying to transform a distance matrix from the qiime2 pathway using a centred log ratio in R. I uploaded the matrix directly from a qiime2 output file, and have converted to a data matrix and have checked that it is symmetrical.

```
structure(list(X1 = c(0, 0.2177609998, 0.2133361674, 0.1549136105, 0.1400395799), X11 = c(0.2177609998, 0, 0.07805820645, 0.1418994689, 0.1934668819), X12 = c(0.2133361674, 0.07805820645, 0, 0.1475390242, 0.1857477705), X13 = c(0.1549136105, 0.1418994689, 0.1475390242, 0, 0.1046740994), X14 = c(0.1400395799, 0.1934668819, 0.1857477705, 0.1046740994, 0)), row.names = c("X1", "X11", "X12", "X13", "X14"), class = "data.frame")
dm <- read.csv(file = "dm.csv", header = TRUE, row.names = 1)
isSymmetric(as.matrix(dm)) [1] TRUE
dm_matrix <- (data.matrix(dm, rownames.force = NA))

```

I then tried to transform the dataset using the `dm_matrix_clr <- data.frame(clr(dm_matrix))` function from the 'compositions' package. This produced an output, however it was not symmetrical. From my understanding, the transformed distance matrix should be symmetric like the input matrix. The same numbers on opposing sides of the diagonal are being transformed, but are resulting in different outputs.

Any help on how resolve this problem so I have a symmetrical transformed distance matrix would be greatly appreciated.

# Answer
Based on conversation in the comments.

**CLR is intended for compositional data. Distance data are not compositional, and compositions are (generally) not symmetric. So you should perform your normalization differently.**

1. If you want to transform your data for normalization purposes on the features themselves (in your case, the organismal relative abundances), you should do it on the compositional input from which you generate your distances. This is the intended application for CLR. For a reference on normalization approaches for compositional data (including CLR), see here.
2. If you want to transform your distance matrix (e.g. the relationships between samples, irrespective of the underlying relative abundances), then you should use one of the recommended methods for transforming distance data.

Note that the applications are very different. Application (1) is transforming the underlying organismal frequency data, such that all of the relative abundances are on the same scale. This will make all downstream analyses more stable / interpretable, including ordination analysis.

Application (2) is saying basically, "my ordination / PCoA on my set of samples *still* looks weird, even after all that normalization I did on my relative abundances. I want to fix it." Examples of "looks weird" given include negative eigenvalues.

You can't substitute (2) for (1). They are working with different data and have different consequences, though both can lead similar issues in your outputs and ordination plots

It seems like your fundamental problem is that qiime2 outputs are difficult to work with. If that is the case, I would suggest asking a separate question about how to handle qiime2 outputs (or contacting the qiime2 authors about how they do compositional normalization).

It looks like the qiime2 forum already has questions about the right way to apply CLR or related transformations in their workflow. Here is another one. I would suggest following those recommendations and spending some quality time with the tutorials and user forum.
> 1 votes

---
Tags: r
---


title: Renaming samples in vcf file

text:
# Question
Title: Renaming samples in vcf file
I am currently asking myself how to rename samples in a vcf file. This is necessary, because the naming of several samples included in thr vcf file was unfortunately not be chosen by me.

I am loading the vcf file via `readVcf()` from the `VariantAnnotation` package in R.

To check the sample names, I can use `samples(header(vcf.file))`, which will give me the names as characters.

It looks like this:

and so on..

(I am actually very sorry that I cannot give a reproducible example here :( )

I also have a vector with all the new names. For example

```
new.names <- c("1","2","3","4")

```

I tried to exchange the "old names" with the new.names vector:

```
samples(header(vcf.file)) <- new.names

```

which leads to the error

```
Error in samples(header(data.vcf)) <- names.sample :
  could not find function "samples<-"

```

I also tried to use the `rutilstimflutre` package, which has a already designed function for this, but unfortunately the package is outdated and not available for the newest R version.

Is there a way to do this? Perhaps a better way to assign the new names to the old ones?

Thanks a lot.

# Answer
Easiest to use `bcftools`, if `samples.txt` has the new sample names:

```
bcftools reheader -s samples.txt -o new.bcf old.bcf

```
> 4 votes

# Answer
This proved to be so simple that it took a while to discover: Just use colnames.

```
x = readVcf("in.vcf")
colnames(x) <- c("new1", "new2")
# At this point, colnames(x) and samples(header(x)) will show different names

# ... but if you write it out, it should use the colnames for both:
writeVcf(x, "out.vcf")
y = readVcf("out.vcf")

# samples(header(y)) and colnames(y) should both be "new1","new2" here

```
> 4 votes

# Answer
Here is an example in python how to rename samples using `pyvcf` submodule I wrote:

```
>>> from fuc import pyvcf
>>> data = {
...     'CHROM': ['chr1', 'chr2'],
...     'POS': [100, 101],
...     'ID': ['.', '.'],
...     'REF': ['G', 'T'],
...     'ALT': ['A', 'C'],
...     'QUAL': ['.', '.'],
...     'FILTER': ['.', '.'],
...     'INFO': ['.', '.'],
...     'FORMAT': ['GT', 'GT'],
...     'Steven': ['0/1', '1/1'],
...     'Sara': ['0/0', '0/1'],
... }
>>> vf = pyvcf.VcfFrame.from_dict([], data)
>>> # vf = pyvcf.VcfFrame.from_file('your_file.vcf')
>>> vf.df
  CHROM  POS ID REF ALT QUAL FILTER INFO FORMAT Steven Sara
0  chr1  100  .   G   A    .      .    .     GT    0/1  0/0
1  chr2  101  .   T   C    .      .    .     GT    1/1  0/1
>>> vf.df = vf.df.rename(columns={'Steven': 'John', 'Sara': 'Rachel'})
>>> vf.df
  CHROM  POS ID REF ALT QUAL FILTER INFO FORMAT John Rachel
0  chr1  100  .   G   A    .      .    .     GT  0/1    0/0
1  chr2  101  .   T   C    .      .    .     GT  1/1    0/1
>>> # vf.to_file('renamed.vcf')

```
> 0 votes

---
Tags: r, vcf
---


# Question
Title: KeyError '/n' in Bio SeqUtils Package ProtParam Script
I am trying to run a function from the ProtParam script in the Bio.SeqUtils package. The function is for calculating the instability index of a protein. The code I am using in my script with the pre-built function is the following:

```
# entry - user inputted string taking in amino acid letters
entry = e1.get()
prot_data = ProteinAnalysis(entry)
result = prot_data.instability_index()

```

When I run the code, I get the following key error related to the instability\_index function.

```
File "/Genome_Toolkit/venv/lib/python3.9/site- packages/Bio/SeqUtils/ProtParam.py", line 152, in instability_index
dipeptide_value = index[this][next]
KeyError: '\n'

```

The instability\_index() function in the ProtParam script looks like this:

```
def instability_index(self):
    index = ProtParamData.DIWV
    score = 0.0

    for i in range(self.length - 1):
        this, next = self.sequence[i: i + 2]
        dipeptide_value = index[this][next]
        score += dipeptide_value

    return (10.0 / self.length) * score

```

I looked at the script myself, specifically the line the error is referring to but couldn't pinpoint the issue. Any guidance is appreciated. Thanks!

the output of `print(entry)` and `print(repr(entry))` looks like this, which is exactly what the entry was in the input:

```
MALWMRLLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGFFYTPKTRREAEDLQVGQVELGG
GPGAGSLQPLALEGSLQKRGIVEQCCTSICSLYQLENYCN
'MALWMRLLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGFFYTPKTRREAEDLQVGQVELGG\nGPGAGSLQPLALEGSLQKRGIVEQCCTSICSLYQLENYCN'

```

# Answer
When the user is pushing `enter` on their keyboard, it is likely inserting a newline character (`\n`) at the end of the user input `entry` and causing the error you have included in your question. You can strip the `\n` by using `entry.rstrip('\n')`:

```
# entry - user inputted string taking in amino acid letters
entry = e1.get()
entry = entry.rstrip('\n')
prot_data = ProteinAnalysis(entry)
result = prot_data.instability_index()

```

---

From the comments where the `print(entry)` was given as

> `MALWMRLLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGFFYTPKTRREAEDLQVGQVELGG GPGAGSLQPLALEGSLQKRGIVEQCCTSICSLYQLENYCN`

and `print(repr(entry))` is given as

> `'MALWMRLLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGFFYTPKTRREAEDLQVGQVELGG\nGPGAGSLQPLALEGSLQKRGIVEQCCTSICSLYQLENYCN'`

the `\n` character is not at the end of the line so `.rstrip('\n')` cannot be used (it only works at the end of the string).

Try instead `.replace('\n', '')`:

```
# entry - user inputted string taking in amino acid letters
entry = e1.get()
entry = entry.replace('\n', '')
prot_data = ProteinAnalysis(entry)
result = prot_data.instability_index()

```

which on my system (removing the `prot_data` and `result` lines) outputs without the `\n` character and `print(repr(entry))` is the same as `print(entry)`.
> 3 votes

---
Tags: proteins, errors, biopython
---


# Question
Title: How to get full alignment output from BLASTP?
I am running a blastp search:

```
time ~/bin/ncbi-blast-2.14.0+-src/c++/ReleaseMT/bin/blastp -db nr -query X.fasta -out nr.X -num_descriptions 4000 -remote

```

but I have 2 problems with the output:

I. Under `Sequences producing significant alignments:`, the output is truncated, and I can't read the exact species name:

```
XP_002906761.1 elongation factor 3, putative [Phytophthora infest...  848     0.0
KAG3108573.1 Elongation factor 3 [Phytophthora idaei] >gb|KAG3151...  845     0.0
KAG2772535.1 Elongation factor 3 [Phytophthora cactorum] >gb|KAG2...  845     0.0
KAG6976665.1 hypothetical protein JG688_00001147 [Phytophthora al...  844     0.0

```

II.

The above sequences do not show up on the detailed sequence alignments below.

I have tried numerous variations of "outfmt" variations, but now am getting warnings that I'm using the system too much. So far, output format `0` is the best.

How can I get the full names and not the truncations?

Why aren't the above species showing up later on, and how can I get the full alignment to show?

# Answer
The problem I had was solved by `num_alignments`, which is more complete than the first section. I do not think that the text truncation is solvable. However, parsing of the latter half solved the problem.

the correct way to get the most complete output is:

`time ~/bin/ncbi-blast-2.14.0+-src/c++/ReleaseMT/bin/blastp -db nr -query X.fasta -out nr.X.0 -outfmt 0 -num_descriptions 4000 -remote -num_alignments 4000`
> -2 votes

# Answer
Hmmm, firstly can you attempt multi-threading on a remote command via  `-num_threads 8`...

```
blastp -db nr -query X.fasta -out nr.X -num_threads 8 -num_descriptions 4000 -remote

```

It's a worth a try. Otherwise it could likely be a function of the network access or remote server resources. I don't think this is a function of the local RAM. Therefore performing a remote blast to download relevant sections of the NCBI database, convert that into a local database and then perform local blast seems the best strategy.

Multi-threading and `-outfmt` work well with local blast.

```
blastp -query query.fa -db test.db -outfmt 6 -max_target_seqs 1000 -num_threads 8 -out blasthits.txt

```

The alignment is a function of the output formats in the above example its `-outfmt 6` that is not useful, but `1` to `4` are relevant as follows.

```
 -outfmt <String>
   alignment view options:
     0 = pairwise,
     1 = query-anchored showing identities,
     2 = query-anchored no identities,
     3 = flat query-anchored, show identities,
     4 = flat query-anchored, no identities,
     5 = XML Blast output,
     6 = tabular,
     7 = tabular with comment lines,
     8 = Text ASN.1,
     9 = Binary ASN.1,
    10 = Comma-separated values,
    11 = BLAST archive format (ASN.1)

```

Multi-threading/parallelisation is the way to overcome memory issues, but really effective only in local operations not remote operations. Again, I think there's too much traffic/ or too many resources requested on NCBI server - and thats where the memory error is occurring.
> 1 votes

---
Tags: blast, blastp
---


# Question
Title: What trajectory analysis method allows to set the form of the trajectory?
I have single cell data from two samples : normal and pathological and I would like to track the progression of the cells (and find the genes that drive it) from normal to pathological.

Upon first glance, it seems most reasonable to define the trajectory **in advance** in such a way that it would have a starting point, an end point, and a straight route between these two. Is it possible, and using what package if so?

P.S. In contrast to the question above - I started with the *destiny* Bioconductor package. However, it always yields three tips for the diffusion pseudotime, and I am not sure how to determine the starting and the end tip of the trajectory...

# Answer
> 4 votes
slingshot (Bioconductor) accepts a priori starts and ends. For differential expression along the trajectory you may check tradeSeq.

---
Tags: bioconductor, single-cell
---


# Question
Title: How to join multiple large .tsv (> 20 GB, bed file) files by two identifier?
*This question has also been asked on Biostars*

I wonder how can I join multiple large .tsv files (over 20 GB each) by using two identifiers (columns) which I used to do using the `full_join` function in `dplyr` package for comparatively smaller files. All files are in `bed format`, where ref\_f.bed will be mapped by the other files by col1 and col2.

```
#merging the files
dflist <- list(ref_f, file_1, file_2, file_3,....file_n)

Test_out <- dflist %>% reduce(full_join, by=c("chr", "start"))

```

```
head ref_f.bed
chr1 1 100 .
chr1 101 200 .
chr1 201 300 .
chr1 301 400 .
chr1 401 500 .
chr1 501 600 .
chr1 601 700 .
chr1 701 800 .
chr1 801 900 .

head file_1.bed
chr1 10101 10200 10158 10177 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10158 10177 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10164 10183 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10164 10183 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10175 10194 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10175 10194 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10175 10194 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10181 10200 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10181 10200 Scerevisiae-UniPROBE-Tbf1.UP00299

head file_2.bed
chr1 10101 10200 10158 10177 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10158 10177 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10164 10183 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10164 10183 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10175 10194 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10175 10194 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10175 10194 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10181 10200 Scerevisiae-UniPROBE-Tbf1.UP00299
chr1 10101 10200 10181 10200 Scerevisiae-UniPROBE-Tbf1.UP00299

head file_3.bed
chr1 29101 29200 29111 29124 Hsapiens-HOCOMOCOv10-AP2D_HUMAN.H10MO.D
chr1 29101 29200 29132 29152 Scerevisiae-UniPROBE-Cha4.UP00300
chr1 29101 29200 29116 29133 Scerevisiae-UniPROBE-Ecm22.UP00305
chr1 29101 29200 29126 29143 Scerevisiae-UniPROBE-Hal9.UP00315
chr1 29101 29200 29126 29143 Scerevisiae-UniPROBE-Hal9.UP00315
chr1 29101 29200 29106 29125 Scerevisiae-UniPROBE-Rpn4.UP00272
chr1 29101 29200 29110 29130 Scerevisiae-UniPROBE-Rsc30.UP00279
chr1 29101 29200 29112 29132 Scerevisiae-UniPROBE-Rsc30.UP00279
chr1 29101 29200 29133 29152 Scerevisiae-UniPROBE-Stp2.UP00338

```

R needs to import all files before performing the task because the file size for more than 20 files surpasses 20 GB, making this task appear difficult for R. This is the first time I have encountered this kind of issue, and I previously used this method in R for merging smaller size files. Therefore, I'm seeking a better solution to this issue.

It would be a great assistance if anyone kindly advised how to overcome this issue.

Thanks

# Answer
> 4 votes
If you are open to a GNU/Linux solution, the following should be a good start:

```
awk '{printf("%s-%s %s %s %s\n", $1, $2, $1, $2, $3, $4);}' ref_f.bed | LC_ALL=C sort -k 1 > sorted_modified_ref_f.txt
awk '{printf("%s-%s %s %s %s %s %s %s\n", $1, $2, $1, $2, $3, $4, $5, $6);}' file_1.bed | LC_ALL=C sort -k 1 > sorted_modified_file_1.txt
join sorted_modified_ref_f.txt sorted_modified_file_1.txt

```

The `awk` commands are creating a "first column" that is a combination of columns 1 and 2 and this is required because `join` "joins" on a single column.

# Answer
> 2 votes
For joining / merging / intersecting bed files, you're better off using bedtools.

It's unclear what you're trying to do with these steps, but I'm guessing it's something like the following:

> The bedtools sub-commands include:
>
> ```
> [ Genome arithmetic ]
>     intersect     Find overlapping intervals in various ways.
>     window        Find overlapping intervals within a window around an > interval.
>     ...
>     merge         Combine overlapping/nearby intervals into a single interval.
>
> ```

If you're interested in combining multiple lines of different files in R, have a look at the `map_dfr` function (reference here), which does a row-wise concatenation of the results of a function mapped to a list.

# Answer
> 0 votes
Thanks for all the responses, however, I got a relatable outcome from (https://unix.stackexchange.com/questions/113898/how-to-merge-two-files-based-on-the-matching-of-two-columns) by using

`$ awk 'NR==FNR {h[$2] = $6; next} {print $1,$2,$3,h[$2]}' file_1 ref_f > out_1`

For `file_1` (`NR==FNR` is only true for the first file argument). Save column 6 in hash-array using column 2 as key: `h[$2] = $6`. Then for `ref_f` and output all three columns `$1,$2,$3`, appending the corresponding saved column from hash-array `h[$2]`.

However, it does not print `col4` and `col5` from `file_1` and must be repeated for each subsequent file. I hope that the answer could be improved further.

---
Tags: r, python, bash, bed, tsv
---


# Question
Title: How can I change all text in a .tsv file to 1, except the first three fields?
I want to convert a .tsv file table in binary format (1 and 0). In which all cells containing text would become `1` and blank cells would become `0`, except first three columns.

The given .tsv file has 6 fields

```
head adrbvbeT_copy.tsv
chr1    1   100
chr1    101 200
chr1    201 300
chr1    301 400
chr1    401 500
chr1    501 600
chr1    601 700
chr1    701 800
chr1    801 900
chr1    901 1000

tail adrbvbeT.tsv
chrY    57226501    57226600
chrY    57226601    57226700    Mmusculus-jaspar2016-Tcf21-MA0832.1
chrY    57226701    57226800    Mmusculus-HOCOMOCOv10-TF7L2_MOUSE.H10MO.A
chrY    57226801    57226900    Mmusculus-HOCOMOCOv10-STAT1_MOUSE.H10MO.A
chrY    57226901    57227000    Mmusculus-HOCOMOCOv10-ZEP2_MOUSE.H10MO.D
chrY    57227001    57227100    NA-SwissRegulon-SMARCC2.SwissRegulon    Mmusculus-HOCOMOCOv10-PBX2_MOUSE.H10MO.C
chrY    57227101    57227200    Mmusculus-HOCOMOCOv10-ZIC3_MOUSE.H10MO.C    Mmusculus-HOCOMOCOv10-TBX5_MOUSE.H10MO.D    Mmusculus-HOCOMOCOv10-PIT1_MOUSE.H10MO.C
chrY    57227201    57227300    Mmusculus-HOCOMOCOv10-MAFB_MOUSE.H10MO.D    Mmusculus-UniPROBE-Tbp.UP00029  Mmusculus-HOCOMOCOv10-MAFB_MOUSE.H10MO.D
chrY    57227301    57227400    Mmusculus-HOCOMOCOv10-PITX2_MOUSE.H10MO.D   Hsapiens-jaspar2016-TGIF2-MA0797.1  Mmusculus-HOCOMOCOv10-PITX2_MOUSE.H10MO.D
chrY    57227401    57227415    Mmusculus-HOCOMOCOv10-TFDP1_MOUSE.H10MO.D   Mmusculus-HOCOMOCOv10-TEAD4_MOUSE.H10MO.D   Mmusculus-HOCOMOCOv10-ZBTB6_MOUSE.H10MO.D

```

First, I converted all blank cells to `0`. Now I'm wondering how I can replace the text in the cells with `1`.

```
awk 'BEGIN { FS = OFS = "\t" } { for(i=1; i<=NF; i++) if($i ~ /^ *$/) $i = 0 }; 1' adrbvbeT_copy.tsv > adrbvbeT_copyN0.tsv

head adrbvbeT_copyN0.tsv
chr1    1   100 0   0   0
chr1    101 200 0   0   0
chr1    201 300 0   0   0
chr1    301 400 0   0   0
chr1    401 500 0   0   0
chr1    501 600 0   0   0
chr1    601 700 0   0   0
chr1    701 800 0   0   0
chr1    801 900 0   0   0
chr1    901 1000    0   0   0

tail adrbvbeT_copyN0.tsv
chrY    57226501    57226600    0   0   0
chrY    57226601    57226700    Mmusculus-jaspar2016-Tcf21-MA0832.1 0   0
chrY    57226701    57226800    Mmusculus-HOCOMOCOv10-TF7L2_MOUSE.H10MO.A   0   0
chrY    57226801    57226900    Mmusculus-HOCOMOCOv10-STAT1_MOUSE.H10MO.A   0   0
chrY    57226901    57227000    Mmusculus-HOCOMOCOv10-ZEP2_MOUSE.H10MO.D    0   0
chrY    57227001    57227100    NA-SwissRegulon-SMARCC2.SwissRegulon    Mmusculus-HOCOMOCOv10-PBX2_MOUSE.H10MO.C    0
chrY    57227101    57227200    Mmusculus-HOCOMOCOv10-ZIC3_MOUSE.H10MO.C    Mmusculus-HOCOMOCOv10-TBX5_MOUSE.H10MO.D    Mmusculus-HOCOMOCOv10-PIT1_MOUSE.H10MO.C
chrY    57227201    57227300    Mmusculus-HOCOMOCOv10-MAFB_MOUSE.H10MO.D    Mmusculus-UniPROBE-Tbp.UP00029  Mmusculus-HOCOMOCOv10-MAFB_MOUSE.H10MO.D
chrY    57227301    57227400    Mmusculus-HOCOMOCOv10-PITX2_MOUSE.H10MO.D   Hsapiens-jaspar2016-TGIF2-MA0797.1  Mmusculus-HOCOMOCOv10-PITX2_MOUSE.H10MO.D
chrY    57227401    57227415    Mmusculus-HOCOMOCOv10-TFDP1_MOUSE.H10MO.D   Mmusculus-HOCOMOCOv10-TEAD4_MOUSE.H10MO.D   Mmusculus-HOCOMOCOv10-ZBTB6_MOUSE.H10MO.D

```

# Answer
> 3 votes
Not an expert, but :

```
awk 'BEGIN { FS = OFS = "\t" } { for(i=4; i<=6; i++) if($i ~ /^ *$/) {$i = 0} else {$i = 1} }; 1'

```

Since you have a very specific file format, this code starts from field 4, checking whether the fields after it are empty or not. I assume also that know that you file has 6 fields max exactly, so I "hardcoded" this knowledge.

---
Tags: bash, bed, awk, text-processing, command-line
---