title: Shell script to validate fastq issue

text:
# Question
Title: Shell script to validate fastq issue
I have two GSM ID as test case where both of them are having data when I checked through sra explorer tool but when I try to fetch through a script only one of them returns output where the other one fails.

**GSM ID** GSM3603268 GSM2683458

**Output** For GSM3603268

```
$ cat GSM3603268.test.fastq
@GSM3603268.1 1 length=100
TTTTTTAAAGAAAACTTGAGCTTTTGGAGCATGGCAACCTAGCCTGCAGACACCGTATCCCCTTGTCCACTTCCCCCTGCAAACCATAAGTCCATTCCTA
+GSM3603268.1 1 length=100
=@@DDDDA++2CDABCEG@GEHHEHGA9EDFHGHD3::B??<DBFFBFB===)5;CA;(7=?;?CE;7)7;>C((583;<9>(,89@@@@CCCCDD:BC:
@GSM3603268.1 1 length=100
AGTTCCACGAGTTTTCTTTTTTTTAAGTGGTAGGAATGGACTTATGGGTTGCAGGGGGAAGTGGACAAGGGGCTACGGTGTCTGCAGGCTTGGTTGCCAT
+GSM3603268.1 1 length=100
:?=BD:D?C)<CDGH4+<CCCGHI)??D<D?00?D*?39C<CG4BFE28;4(.7AB############################################

$ cat GSM3603268.test.log
2023-09-26T21:18:46 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
2023-09-26T21:18:46 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
2023-09-26T21:18:46 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
Read 1 spots for GSM3603268
Written 1 spots for GSM3603268

```

Output for GSM2683458

```
$ cat GSM2683458.test.fastq
    bash-5.1#

$ cat GSM2683458.test.log
2023-09-26T21:16:55 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
2023-09-26T21:16:55 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
2023-09-26T21:16:55 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
2023-09-26T21:17:11 fastq-dump.3.0.1 warn: directory not found while opening manager within virtual file system module - can't open NC_000076.5 as a RefSeq or as a WGS
2023-09-26T21:17:11 fastq-dump.3.0.1 err: directory not found while opening manager within virtual file system module - failed GSM2683458

=============================================================
An error occurred during processing.
A report was generated into the file '/root/ncbi_error_report.txt'.
If the problem persists, you may consider sending the file
to 'sra-tools@ncbi.nlm.nih.gov' for assistance.
=============================================================

fastq-dump quit with error code 3

```

The script

```
#!/bin/sh

set -x
PS4='[\\d \\t] '

# Check parameter for error
check=0
# Print fastq-dump executable path
echo $(which fastq-dump)

# Function to download FastQ with retries
download_with_retry() {
    local id="$1"
local retries=5
for attempt in $(seq "$retries"); do
    # Download start of fastq
    fastq-dump $(get_ngc) -X 1 -Z --split-spot "$id" > "${id}.test.fastq" 2> "${id}.test.log"
    numLines=$(cat "${id}.test.fastq" | wc -l)
    if [ $numLines -gt 0 ]; then
            echo "${id} has data... OK"
        return 0
    else
        echo "${id} does not have data on attempt $attempt... Retrying in 5 seconds..."
        sleep 5
    fi
done
echo "${id} could not be downloaded after $retries attempts... ERROR"
    check=1
    return 1
}

# Loop through all parameters to check validity
for file in "$@"; do
cp "${file}" .
    # Extract filename for sampleID
    file_basename=$(basename "${file}")
    id="${file_basename%".id"}"
# Start validation with retries
echo "Checking ${id}..."
    download_with_retry "$id"
done

# Exit with error if some FastQs are not accessible
if [ $check -gt 0 ]; then
    echo "ERROR: One or more samples have inaccessible FastQs.. exiting"
    exit 1
fi

```

What exactly is going wrong in the case of the other ID? No matter how many retries I give, it fails.

Any suggestion or help would be really appreciated.

The error log:

```
$ cat ncbi_error_report.txt
<Report>
 <Run>
  <Date>
   <Start value="Tue Sep 26 2023 9:16:55 PM"/>
   <End value="Tue Sep 26 2023 9:17:11 PM"/>
  </Date>
  <Home name="HOME" value="/root"/>
  <Cwd>/tmp</Cwd>
  <CommandLine argc="6">
   <Arg index="0" value="fastq-dump"/>
   <Arg index="1" value="-X"/>
   <Arg index="2" value="1"/>
   <Arg index="3" value="-Z"/>
   <Arg index="4" value="--split-spot"/>
   <Arg index="5" value="GSM2683458"/>
  </CommandLine>
  <Result rc="RC(rcVFS,rcMgr,rcOpening,rcDirectory,rcNotFound)"/>
  <User admin="true"/>
 </Run>
 <Configuration>
  <Files count="2">
   <File name="/etc/ncbi/settings.kfg"/>
   <File name="/root/.ncbi/user-settings.mkfg"/>
  </Files>
  <refseq state="not found"/>
  <krypto state="pwfile: not found"/>
  <sra>
   <quality_type>raw_scores</quality_type>
  </sra>
  <Config>
  <ConfigurationFiles>
/etc/ncbi/settings.kfg
/root/.ncbi/user-settings.mkfg
    </ConfigurationFiles>
    <APPNAME>"fastq-dump"</APPNAME>
    <APPPATH>"/tmp/"</APPPATH>
    <BUILD>"RELEASE"</BUILD>
    <HOME>"/root"</HOME>
    <HOST></HOST>
    <LIBS>
      <GUID>"119c217a-7b81-47e8-91d6-56d19c8c9f15"</GUID>
      <IMAGE_GUID>"119c217a-7b81-47e8-91d6-62229c64ee59"</IMAGE_GUID>
    </LIBS>
    <NCBI_HOME>"/root/.ncbi"</NCBI_HOME>
    <NCBI_SETTINGS>"/root/.ncbi/user-settings.mkfg"</NCBI_SETTINGS>
    <OS>"linux"</OS>
    <PWD>"/tmp"</PWD>
    <USER></USER>
    <VDB_CONFIG></VDB_CONFIG>
    <VDB_ROOT></VDB_ROOT>
    <kfg>
      <arch>
        <bits>"64"</bits>
        <name>"56d19c8c9f15"</name>
      </arch>
      <dir>"/root/.ncbi"</dir>
      <name>"user-settings.mkfg"</name>
    </kfg>
    <libs>
      <cloud>
        <report_instance_identity>"true"</report_instance_identity>
      </cloud>
    </libs>
    <repository>
      <user>
        <ad>
          <public>
            <apps>
              <file>
                <volumes>
                  <flat></flat>
                  <flatAd>"."</flatAd>
                </volumes>
              </file>
              <refseq>
                <volumes>
                  <refseqAd>"."</refseqAd>
                </volumes>
              </refseq>
              <sra>
                <volumes>
                  <sraAd>"."</sraAd>
                </volumes>
              </sra>
              <sraPileup>
                <volumes>
                  <ad>"."</ad>
                </volumes>
              </sraPileup>
              <sraRealign>
                <volumes>
                  <ad>"."</ad>
                </volumes>
              </sraRealign>
              <wgs>
                <volumes>
                  <wgsAd>"."</wgsAd>
                </volumes>
              </wgs>
            </apps>
            <root>"."</root>
          </public>
        </ad>
      </user>
    </repository>
    <sra>
      <quality_type>"raw_scores"</quality_type>
    </sra>
    <vdb>
      <lib>
        <paths>
          <kfg>"/usr/local/bin"</kfg>
        </paths>
      </lib>
    </vdb>
  </Config>
  <RemoteAccess available="false"/>
  <CurrentProtectedRepository found="false"/>
 </Configuration>
 <Object path="https://sra-pub-run-odp.s3.amazonaws.com/sra/SRR5755657/SRR5755657" type="database" fs_type="unexpected">
  <Dependencies>
   <List count="22" missing="22">
    <Dependency index="0" seq_id="NC_000067.5" local="false" path=""/>
    <Dependency index="1" seq_id="NC_000068.6" local="false" path=""/>
    <Dependency index="2" seq_id="NC_000069.5" local="false" path=""/>
    <Dependency index="3" seq_id="NC_000070.5" local="false" path=""/>
    <Dependency index="4" seq_id="NC_000071.5" local="false" path=""/>
    <Dependency index="5" seq_id="NC_000072.5" local="false" path=""/>
    <Dependency index="6" seq_id="NC_000073.5" local="false" path=""/>
    <Dependency index="7" seq_id="NC_000074.5" local="false" path=""/>
    <Dependency index="8" seq_id="NC_000075.5" local="false" path=""/>
    <Dependency index="9" seq_id="NC_000076.5" local="false" path=""/>
    <Dependency index="10" seq_id="NC_000077.5" local="false" path=""/>
    <Dependency index="11" seq_id="NC_000078.5" local="false" path=""/>
    <Dependency index="12" seq_id="NC_000079.5" local="false" path=""/>
    <Dependency index="13" seq_id="NC_000080.5" local="false" path=""/>
    <Dependency index="14" seq_id="NC_000081.5" local="false" path=""/>
    <Dependency index="15" seq_id="NC_000082.5" local="false" path=""/>
    <Dependency index="16" seq_id="NC_000083.5" local="false" path=""/>
    <Dependency index="17" seq_id="NC_000084.5" local="false" path=""/>
    <Dependency index="18" seq_id="NC_000085.5" local="false" path=""/>
    <Dependency index="19" seq_id="NC_000086.6" local="false" path=""/>
    <Dependency index="20" seq_id="NC_000087.6" local="false" path=""/>
    <Dependency index="21" seq_id="NC_005089.1" local="false" path=""/>
   </List>
  </Dependencies>
 </Object>
 <SOFTWARE>
  <VDBLibrary vers="2.7.47"/>
  <Build static="true">
   <Module name=""/>
  </Build>
  <Tool date="Nov 18 2022" name="fastq-dump" vers="3.0.1">
   <Binary path="/usr/local/bin/fastq-dump" type="alias" md5="c461c39bfa514aff3c4f7c0416ced617">
    <Alias resolved="fastq-dump.3">
     <Alias resolved="fastq-dump.3.0.1">
      <Alias resolved="sratools.3.0.1"/>
     </Alias>
    </Alias>
   </Binary>
  </Tool>
 </SOFTWARE>
 <Env>
 </Env>
</Report>

```

# Answer
Use the SRR accessions, not the sample-IDs you are using.

With entrez e-utils, I found the SRR of GSM2683458:

`esearch -db sra -query GSM2683458 | efetch -format xml | xtract -pattern IDENTIFIERS -if PRIMARY_ID -contains SRR -element PRIMARY_ID` is `SRR5755657`.

For GSM3603268, using `esearch -db sra -query GSM3603268 | efetch -format xml | xtract -pattern IDENTIFIERS -if PRIMARY_ID -contains SRR -element PRIMARY_ID`, the SRR is: `SRR8571942`

Then just do `fasterq-dump SRR5755657; fasterq-dump SRR8571942`.
> 2 votes

---
Tags: fastq, sratoolkit, ncbi
---


# Question
Title: Compare files with similar names using placeholders/wildcards in Snakemake
In Snakemake I want to use two input files for comparison. The one file has a sample name (e.g. 526846), with an experimental id (e.g. EF\_S13' and 'MF\_S5'). The file I want to compare it to has the same sample name only. I want to compare each sample with experiment ID (first files) to the file with corresponding sample name (second file). These files are in different directories. I thought I would be able to use two expand functions with two lists, as shown in the code below. Note: these are example lists and the lists in reality is muchlarger.

Relevant code:

```

SAMPLE_SET = ['12345', '56789', '54321']
ID_SET = ['EF_S13', 'MF_S5', 'EF_S11', 'MF_S3', 'EF_S10', 'MF_S2']

rule all:
    input:
        array = expand("results/{sample}_{samp_id}_array/{sample}_{samp_id}.summary.csv", zip, sample = SAMPLE_SET, samp_id = ID_SET)

rule happy_rtg_array:
    """ Normalizes and compares vcf files with array samples"""
    input:
        query_vcf = expand("WESP_run/results/variant_calling/{sample}_{samp_id}.vcf", zip, sample = SAMPLE_SET, samp_id = ID_SET),
        truth_vcf = expand(os.path.join(config["truth_dir_vcf_gsa"],"{sample}.vcf.gz"), sample = SAMPLE_SET),
    output:
        "WESPValidation/results/{sample}_{samp_id}_array/{sample}_{samp_id}.summary.csv"
[more code]

```

However, when I look at the output log file with --dry-run, I get as input for the job a huge list with all the inputs in one go, rather then per file/comparison, and one output file (see relevant output), and that times the amount of the (right) combinations of the file name and id (counted from the output).

Also the wildcards are not what I thought they should be. I might be confused with the expand function, as I thought how I used them now should "lock in" the wildcards and become placeholders filled in by the lists.

Relevant output:

```
[Tue Oct  3 14:31:20 2023]
rule happy_rtg_array:
    input: 12345_EF_S13.vcf, 12345_MF_S5.vcf, 56789_EF_S11.vcf, 56789_MF_S3.vcf, 54321_EF_S10.vcf, 54321_MF_S2.vcf, 12345.vcf, 56789.vcf, 54321.vcf
    output: results/12345_EF_S13_array/12345_EF_S13.summary.csv
    log: happy_rtg_array/12345_EF_S13.log
    jobid: 6
    reason: Missing output files: WESPValidation/results/12345_EF_S13_array/12345_EF_S13.summary.csv
    wildcards: sample=12345_EF, samp_id=S13
    threads: 10
    resources: mem_mb=10000, disk_mb=7855, tmpdir=/tmp, threads=1, runtime_min=120

```

I doubt this is how it should be, and wonder what I am doing wrong. Is there a better/correct way to do this.

Thank you in advance

# Answer
After some trial and error, I fixed it with the expand function only in the rule all and removing the expand from the input in the rule happy\_rtg\_array (thanks to @Pallie), and by changing the output file paths:

```
SAMPLE_SET = ['12345', '56789', '54321']
ID_SET = ['EF_S13', 'MF_S5', 'EF_S11', 'MF_S3', 'EF_S10', 'MF_S2']

rule all:
    input:
        array = expand("results/{sample}/{samp_id}/summary.csv", zip, sample = SAMPLE_SET, samp_id = ID_SET)

rule happy_rtg_array:
    """ Normalizes and compares vcf files with array samples"""
    input:
        query_vcf = "_".join(["WESP_run/results/variant_calling/{sample}","{samp_id}.vcf"]),
        truth_vcf = os.path.join(config["truth_dir_vcf_gsa"],"{sample}.vcf.gz"),
    output:
        "WESPValidation/results/{sample}/{samp_id}/summary.csv"
[more code]

```

Basically, I replaced all the underscores with slashes. For some reason these underscores are messing up the wildcards when trying to connect them with it. Maybe someone more knowledgeable knows why this is.

```
[Tue Oct  5 12:24:10 2023]
rule happy_rtg_array:
    input: 12345_EF_S13.vcf
    output: results/12345/EF_S13/summary.csv
    log: happy_rtg_array/12345_EF_S13.log
    jobid: 6
    reason: Missing output files: WESPValidation/results/12345/EF_S13/summary.csv
    wildcards: sample=12345, samp_id=EF_S13
    threads: 10
    resources: mem_mb=10000, disk_mb=7855, tmpdir=/tmp, threads=1, runtime_min=120

```
> 0 votes

# Answer
Using `expand` in your rule's `input:` causes snakemake to aggregate all vcfs into one list and leads to your undesired behaviour. By simply removing the `expand` calls from the `input:` you should get your desired result.

https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#the-expand-function
> 1 votes

---
Tags: python, snakemake
---


title: How to interpret a GenomeScope plot?

text:
# Question
Title: How to interpret a GenomeScope plot?
I am working on a new sponge genome and I have produced the genomescope plot from the 21-mer kmer frequencies. I have hard time interpreting the plot. Can someone please help me? Thank you.

# Answer
Based on the numbers shown above the plot, and the appearance of the curve, this looks like reads from a diploid genome approximately 113Mb in length that was sequenced at around 174X coverage.
> 1 votes

---
Tags: ngs, genome, sequence-annotation, illumina, genome-sequencing
---


# Question
Title: How to create Phylogenetic Trees from fasta files in Python or R?
I have around a hundred Fasta files (and will collect several thousand) with DNA sequences and +50x coverage. What is a recommended method to construct a phylogenetic tree? Solutions in Python or R are sought.

I found `Phylo` from `Biopython` only handles already calculated trees.

# Answer
I would not look for a package for this, but instead build a small pipeline calling external tools with something like the following workflow:

* Cluster the ~100 sequences with CD-HIT-EST/PSI-CD-HIT or many other options
* Take all the sequences that form one individual cluster and build a multiple sequence alignment (MSA) with MAFFT/ClustalOmega or similar
* Take the MSA and build a phylogenetic tree with a Maximum-Likelihood approach like IQ-TREE or similar
* Visualize the tree file with Jalview or similar

Of course this is rather general and depending on exactly what you're doing you may want a different workflow and/or different tools. You should also explore the parameter space, *do not* assume the defaults are necessarily good choices
> 7 votes

# Answer
The obvious single answer is R "ape". This will give you access to PhylML for tree building and Clustal/Muscle for alignment building. The paths to the binarys are important. There are several distance methods in there such as NJ and BIONJ. Its distance approaches however don't look mainstream, but I could be wrong.

There are functions within ape which are cool, the tree sorting is very cool and I need to read through this with much greater care. Personally I wouldn't perform a core phylogenetic analysis within R, because the standalones are sufficient and the analysis is intensive.

https://cran.r-project.org/web/packages/ape/ape.pdf
> 3 votes

# Answer
I agree with Chris Rands that a reasonable approach would be to call external tools.

However, if you really want to do the phylogeny from within Python, you could use the P4 package, which is a bit complicated to handle but gives you lots of options in the way to build MCMC-based bayesian phylogenies:

https://github.com/pgfoster/p4-phylogenetics

You would still need something else to align the sequences before.

To visualize the tree using python, you could use the ete toolkit, which is likely more powerful than what you can find in Biopython: http://etetoolkit.org/
> 3 votes

# Answer
(if I understand your situation correctly) `https://www.rdocumentation.org/packages/seqinr/versions/3.6-1/topics/read.alignment` shows how to use the function **read.alignment** which can take fasta msf etc. The docs provide the example' `read.alignment(file = system.file("sequences/LTPs128_SSU_aligned_First_Two.fasta", package = "seqinr"), format = "fasta", whole.header = TRUE)` but you can use this code below (assumes those files are aligned) to go from reading the tree to getting the distances, producing the neighbor joining phylogenetic tree, and then plotting the tree.

```
library("Biostrings")
library("seqinr")
library("ape")
library(phylogram)
library("dendextend")

fasta.res <- read.alignment(file = "geneticAlignment.msf", format = "fasta")
fasta.res.dist.alignment = dist.alignment(msf.res, matrix = "identity")
fasta.res.dist.alignment.nj = nj(fasta.res.dist.alignment)
plot(fasta.res.dist.alignment.nj, main = "from fasta files")

```
> 2 votes

# Answer
Within the R/Bioconductor ecosystem you could use the `DECIPHER` package (https://bioconductor.org/packages/release/bioc/html/DECIPHER.html) to align your sequences and then continue with either `ape` or `phangorn` (https://klausvigo.github.io/phangorn/) packages to compute a tree. You can then visualize your tree using for example the `ggtree` package (https://bioconductor.org/packages/release/bioc/html/ggtree.html).
> 1 votes

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


# Question
Title: Why do we delete scaffolds shorter than 500 - 1000 bp from the assembled genome?
After assembly of genome, some protocols sometimes call for removal of scaffolds shorter than 500 or 1000 (some papers have one number while the other has the other.) Is this simply to remove the artefact and to make later analysis smoother, without the extroneous data? Also how do we know how where is our cut off point for size?

# Answer
I recognise my answer is late but I hope it helps. The following has no emperical proof, just anecdotal from working with bacterial genomes:

I'm not sure what kind of genome you're working with, but from my work with bacteria, we often delete scaffolds shorter than 1000 - 3000 bp because these are considered "noise". Essentially the scaffolds are too small to gather any particularly useful information from as they may only carry a single gene (or a handful of little genes). If you've added in long read data with your short, half the time those small scaffolds are just artifacts from sequencing errors.
> 2 votes

---
Tags: assembly
---


# Question
Title: How can I retrieve the GRCh38 coordinates of a list of rsids?
I have a list of about 100,000 rsids and I want to get their genomic coordinates on the GRCh38 genome build. Is there a command line tool that allows me to do this? If yes, which one?

I have tried bcftools but, given the error message I got, I believe that it does not support search by rsid:

```
bcftools view -R rsids_chr1_subset.txt -f PASS -v snps dbSNP_v156_GRCh38.p14.vcf.gz > rsids_chr1.vcf
[E::bcf_sr_regions_init] Could not parse 1-th line of file rsids_chr1_subset.txt, using the columns 1,2[,-1]
Failed to read the regions: rsids_chr1_subset.txt

```

# Answer
Several years ago I created rsidx to address this particular use case. https://github.com/bioforensics/rsidx

```
# VCF should be sorted by genomic coordinates and indexed by tabix
rsidx index dbSNP151_GRCh38.vcf.gz dbSNP151_GRCh38.rsidx
rsidx search dbSNP151_GRCh38.vcf.gz dbSNP151_GRCh38.rsidx rs3114908 rs10756819

```
> 5 votes

# Answer
rsids\_chr1.txt contains a list of rsids, one per line.

```
bcftools view --include ID==@rsids_chr1.txt dbSNP_v156_GRCh38.p14.vcf.gz -Ou | bcftools query -f '%CHROM  %ID %POS  %REF  %ALT{0}\n' > dbSNP_v156_GRCh38.p14_AF.vcf

```
> 2 votes

# Answer
This was part of a custom liftover function which handled multiple alleles and reverse compliments. I do not know if it will be useful. I believe there was a separate thing somewhere with lookup that I may add. It handled apparent non-1:1 mappings in addition to platform-specific "rsIDs."

converter = get\_lifter(current, target)

```
def convert_coordinates(input_file, output_file, converter):
        line_count = 0
        skipped = 0
        strand_flips = 0
        with open(input_file, 'r') as infile, open(output_file, 'w') as outfile:
            for line in infile:
                # Skip headers
                if line.startswith("#"):
                    outfile.write(line)
                    continue

                # VCF files are tab-separated
                parts = line.strip().split('\t')

                # Get the chromosome and position
                chrom = parts[0]
                pos = int(parts[1])

                # Convert the coordinates
                try:
                    new_coords = converter.query(chrom, pos)
                except KeyError:
                    #print(f"Chromosome {chrom} not found in the liftover chain file.")
                    skipped += 1
                    continue

                # If there's no new coordinates, we skip this variant
                if not new_coords:
                    skipped += 1
                    continue

                # Get the new chromosome, position and strand
                new_chrom, new_pos, strand = new_coords[0]

                # If the new strand is negative, reverse complement the REF and ALT alleles
                if strand == '-':
                    parts[3] = reverse_complement(parts[3])
                    parts[4] = reverse_complement(parts[4])
                    strand_flips += 1

                # Replace the old chromosome and position with the new ones
                parts[0] = new_chrom
                parts[1] = str(new_pos)

                # Write the line back into the file
                outfile.write('\t'.join(parts))
                outfile.write('\n')

                # Print a progress statement every 1000 lines
                line_count += 1
                if line_count % 100000 == 0:
                    print(f"Processed {line_count} lines")

        # Calculate and print the percentage of variants skipped
        percentage_skipped = (skipped / line_count) * 100
        percentage_flips = (strand_flips / line_count) * 100
        print(f"\nFinished processing. Out of {line_count} total lines:")
        print(f"{percentage_skipped}% of variants were skipped due to missing chromosomal data.")
        print(f"{percentage_flips}% of variants required a strand flip.")

```
> 1 votes

---
Tags: bcftools, rsid
---


# Question
Title: Can the PGS Calculator give the contribution of each variant to the overall PGS?
I am using the PGS Calculator (https://github.com/PGScatalog/pgsc\_calc/) to compute the polygenic scores of individuals based on a scoring file.

PGS Calculator works wonderfully and outputs the aggregated score for each individual.

However, I am interested in the contribution of each variant to the score of each individual.

At the moment, I have:

| IID | PGS |
| --- | --- |
| 33432984732 | 10 |
| 32432809751 | 20 |

I would like the following table:

| IID | Variant 1 | Variant 2 | Variant 3 |
| --- | --- | --- | --- |
| 33432984732 | 3 | 4 | 3 |
| 32432809751 | 6 | 8 | 6 |

in which I can see the contribution of each variant to the overall score.

Is this possible with PGS Calculator? If not, is there a tool other than PGS Calculator that I can use? Or should I use my own scripts?

# Answer
> 0 votes
This would be a hard feature to support natively, because some PRSs have 10 million+ SNPs.

There is a tool that does this, but I don't want to recommend it because it will require debugging their source code upon running, and it is less accurate than what you are using here.

I suggest doing --scorefile "my\_custom\_scores/\*.txt" if you are interested in only a few variants. Add only one score to this file. This may cause some error especially if some of your samples do not have all positions defined.

If you have more variants, there is a way to combine multiple scores in a single scorefile instead of a directory.

---
Tags: polygenic-risk-score
---


# Question
Title: How to fix "No gene can be mapped, Expected input gene ID, return NULL..." error in clusterProfiler?
I am trying EnrichKEGG on a few gene samples of genes to check it works before inputting the full genome:

This is my data input as test.csv

```
id,fc
AJAP_14870,-0.04
AJAP_14875,0.32
AJAP_14880,-0.06
AJAP_14885,0.15
AJAP_14890,-0.52
AJAP_14895,-0.39
AJAP_14900,0.16

```

```
geneList = test[,2]
names(geneList) = as.character(test[,1])
geneList = sort(geneList, decreasing = TRUE)

mygene <- names(geneList)[abs(geneList) > 2]

kk <- enrichKEGG(gene         = mygene,
                 organism     = 'aja',
                 pvalueCutoff = 0.05)

```

There are no spaces in the gene names, I tested via

```
print(head(test$ID), quote=TRUE)

```

But I get the error

```
--> No gene can be mapped....
--> Expected input gene ID:
--> return NULL...

```

I also tried installing through other means

```
getOption("clusterProfiler.download.method")
#[1] "libcurl"
install.packages("R.utils")
library(R.utils)
R.utils::setOption("clusterProfiler.download.method","auto")

```

but that gives me the same error when I run the code. Any ideas how to fix this?

Thanks

# Answer
> 3 votes
This issue may arise with earlier versions of clusterProfiler.

A solution is to re-install clusterProfiler : `BiocManager::install("clusterProfiler", force = TRUE)`

Pay attention that updating your R or Bioconductor (and restart your R session) version may be needed as noted in the comments

# Answer
> 1 votes
`enrichKEGG()` expects entrez gene IDs so you will need to convert your gene identifiers to entrez identifiers first.

EDIT: Please see the answer by @Basti.

---
Tags: r, bioconductor, clusterprofiler
---


title: Extract data from genome

text:
# Question
Title: Extract data from genome
*This question was also asked on Biostars 3 months ago, i.e. a long time - so its perfectly reasonable for the question to be asked here*

I started studying bioinformatics and I enjoy it a lot (I'm a software developer with around 15 years of experience but mostly in backend, scripting, architecture, and security) and I want to do, as a project, a database regarding natural compounds. When I saw that on NCBI there are a lot of plant genomes I thought, that, with software, somehow we can "guess" (probabilities, not certainties) what a plant can contain or not.

I'm wondering, if, by bioinformatic means only, it is possible to "extract" genes and, after that proteins and natural compounds from plant genome data that is available on NCBI (please see https://www.ncbi.nlm.nih.gov/datasets/genome/GCA\_029618835.1/ as an example)? Example: just by processing the data from the link above we can get to Flavoxanthin or any other substance.

I was wondering if a kind of "de-novo" annotation can be made (if that's the correct term).

What other useful information from a chemical compound perspective can be found from digging on the genome itself?

Thank you!

# Answer
This sounds like you are wanting to perform a machine learning calculation. However, for this question leveraging KEGG, Kyoto Encyclopaedia of Genome and Gene and building an SQL database around it is much more effective. It represents around 60% of all functional genes. KEGG is the best gene annotation out there IMO, particularly for flavonoids etc... Fantastic and all their linking pathways, everything links to everything, just amazing. KEGG represents an enormous volume of work, much of it being manually crunching publications for each gene.

The entire database of KEGG keys is downloadable. Connecting a gene to a class of function is easy via KEGG the relationships are defined its just SQL-ing the whole thing together (querying all the tables). Again, it's totally and completely brilliant. Sure they'll be stuff out there has SQL-ed KEGG key tables, but if your a developer its better to build it independently - this way you're certain whats in there and what its doing.

However, linking raw sequence data to a KEGG gene code is more challenging. lastKOALA and GhostKOALA is ... yeah well an approach no more. There's a main stream approach but with current limitations in implementation. KEGG however isn't limiting, it's linking a gene to a KEGG id is the challenge. Getting hold of the raw sequence for KEGG would be great, fantastic in fact, but I'm yet to find it.
> 2 votes

# Answer
I've used Augustus stand alone (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1160219/) to find genes in a given eukaryotic genome. I agree you can then use things like KofamSCAN to categorize the genes against the KEGG database. If you want to compare against the NCBI database, diamond blast (https://github.com/bbuchfink/diamond) is probably the quickest way to go :)
> 1 votes

---
Tags: genome, gene, sequence-annotation, genomics, plants
---


title: Compare substitution matrices

text:
# Question
Title: Compare substitution matrices
I am looking for methods to compare substitution matrices to determine which one is better for a given set of data / type of protein.
Do you know such methods / tools?

e.g. BLOSUM62 vs PHAT for 50k - 100k+ sequences

# Answer
The technique is phylogenetics is called likelihood ratio and for a regular amino acid matrix is 2Δ*ln* which is assessed on Chi2 tables where Δ*ln* is difference in likelihood between model 1 and model 2. The degrees of freedom hmmmm... is the number of parameters that are constrained, this involves the number of taxa minus the number of parameter differences between the models. This makes performing 50 000 taxa ... tricky because the degrees of freedom will be outside the Chi2 table.

iqtree will perform a tree and provide a likelihood for 50+K taxa. There are two limitations,

* No PHAT model
* degrees of freedom (as mentioned)

Overcoming df is simply to say whatever is the lowest likelihood is the best model and simply state Δ*ln* without bothering with significance testing. If the difference between models is large you'll get away with that (needs to be big) and for 50+K taxa it will be need to be huge.

Getting PHAT (the model below the ---- might be PHAT) into iqtree is certainly doable because of QMaker. However, how this is done I don't know, because I never done this calculation with iqtree, you'll need to read the manual. For traditional packages such as PAML it's easy ... you make a tree then specific the custom amino acid model and it spits out the likelihood. For the number of taxa involved, PAML is not possible, it's impossibly big for PAML, even 1000 is a lot of PAML.

iqtree does have LG, WAG and JTT amino acid models and will almost certainly have all the BLOSUM matrices and you just specify those at the command line. Again a comparison with PHAT needs to specify a custom matrix and I don't know how to do that in iqtree. The only thing I know is it must be doable.

iqtree is the most efficient of the maximum likelihood phylogenetic methods by a very long way. Hence for dealing with 50+K is why I'm recommending it here.

**Note** The theory about whats going on and why the whole thing needs phylogenetics is very complicated. Basically its an entire lecture on phylogenetic theory.

---

This I believe is PHAT (needs checking)

https://github.com/biopython/biopython/pull/283

Note the code is copied and pasted from the above link, therefore it should be strictly regarded as a quote for copyright purposes.

```
#!/usr/bin/env python
matrix = {
('W', 'F') :  0.0, ('L', 'R') : -6.0, ('I', 'I') :  5.0, ('Q', 'Q') :  9.0,
('W', 'N') : -5.0, ('V', 'I') :  3.0, ('H', 'T') : -4.0, ('H', 'P') : -6.0,
('W', 'V') : -4.0, ('Q', 'E') :  1.0, ('W', 'R') : -7.0, ('Q', 'A') : -3.0,
('H', 'H') : 11.0, ('H', 'D') : -1.0, ('L', 'N') : -3.0, ('Y', 'M') : -2.0,
('Y', 'I') : -3.0, ('Y', 'E') : -2.0, ('E', 'S') : -3.0, ('Y', 'A') : -3.0,
('Y', 'Y') : 11.0, ('T', 'C') : -1.0, ('E', 'C') : -7.0, ('Y', 'Q') :  0.0,
('E', 'G') : -3.0, ('V', 'A') :  1.0, ('C', 'C') :  7.0, ('M', 'R') : -6.0,
('P', 'T') : -4.0, ('V', 'E') : -5.0, ('P', 'P') : 13.0, ('I', 'T') : -1.0,
('K', 'S') : -5.0, ('R', 'G') : -5.0, ('I', 'P') : -4.0, ('R', 'C') : -8.0,
('A', 'T') :  0.0, ('K', 'K') :  5.0, ('A', 'P') :  3.0, ('V', 'M') :  1.0,
('I', 'D') : -5.0, ('K', 'C') :-10.0, ('K', 'G') : -5.0, ('R', 'S') : -6.0,
('F', 'Q') : -2.0, ('F', 'A') : -1.0, ('V', 'V') :  4.0, ('M', 'N') : -2.0,
('F', 'E') : -5.0, ('D', 'N') :  2.0, ('F', 'I') :  0.0, ('F', 'M') :  0.0,
('M', 'S') : -2.0, ('S', 'S') :  6.0, ('L', 'Q') : -3.0, ('W', 'E') : -7.0,
('W', 'A') : -4.0, ('W', 'M') : -4.0, ('H', 'S') : -2.0, ('W', 'I') : -4.0,
('S', 'C') :  1.0, ('L', 'A') : -1.0, ('L', 'E') : -5.0, ('W', 'Q') :  1.0,
('H', 'G') : -4.0, ('Q', 'N') :  2.0, ('H', 'C') : -7.0, ('L', 'M') :  2.0,
('W', 'Y') :  1.0, ('Y', 'N') :  2.0, ('E', 'P') : -5.0, ('Y', 'F') :  4.0,
('E', 'T') : -5.0, ('A', 'A') :  5.0, ('I', 'N') : -3.0, ('G', 'A') :  1.0,
('Y', 'V') : -3.0, ('E', 'D') :  6.0, ('W', 'H') : -3.0, ('Y', 'R') : -6.0,
('M', 'Q') : -1.0, ('P', 'S') : -3.0, ('R', 'H') : -4.0, ('A', 'C') :  1.0,
('R', 'D') : -7.0, ('K', 'P') : -4.0, ('L', 'D') : -5.0, ('K', 'T') : -6.0,
('V', 'N') : -3.0, ('M', 'A') : -1.0, ('K', 'H') : -5.0, ('V', 'R') : -7.0,
('P', 'C') : -8.0, ('M', 'E') : -5.0, ('A', 'S') :  2.0, ('T', 'T') :  3.0,
('R', 'T') : -6.0, ('I', 'G') : -2.0, ('R', 'P') : -7.0, ('K', 'D') : -5.0,
('I', 'C') : -3.0, ('F', 'R') : -7.0, ('F', 'V') : -1.0, ('L', 'C') : -2.0,
('F', 'F') :  6.0, ('D', 'A') : -5.0, ('F', 'N') : -1.0, ('W', 'D') : -7.0,
('L', 'P') : -5.0, ('Q', 'S') : -1.0, ('N', 'C') : -2.0, ('N', 'G') : -1.0,
('H', 'N') :  4.0, ('W', 'T') : -7.0, ('Q', 'G') : -2.0, ('W', 'P') : -6.0,
('Q', 'C') : -5.0, ('N', 'S') :  1.0, ('L', 'H') : -4.0, ('L', 'L') :  4.0,
('G', 'T') : -1.0, ('M', 'M') :  6.0, ('G', 'P') : -3.0, ('Y', 'K') : -4.0,
('Y', 'G') : -3.0, ('Y', 'C') : -1.0, ('E', 'A') : -5.0, ('E', 'E') : 12.0,
('Y', 'S') : -2.0, ('M', 'P') : -5.0, ('V', 'C') : -2.0, ('M', 'T') :  0.0,
('V', 'G') : -2.0, ('R', 'E') : -6.0, ('V', 'K') : -8.0, ('K', 'Q') : -1.0,
('R', 'A') : -6.0, ('I', 'R') : -6.0, ('N', 'A') : -2.0, ('V', 'S') : -2.0,
('M', 'D') : -5.0, ('M', 'H') : -4.0, ('K', 'A') : -7.0, ('R', 'Q') : -2.0,
('K', 'E') : -4.0, ('F', 'S') : -2.0, ('I', 'K') : -7.0, ('D', 'P') : -5.0,
('D', 'T') : -5.0, ('I', 'M') :  3.0, ('F', 'C') :  0.0, ('W', 'L') : -3.0,
('F', 'G') : -2.0, ('F', 'K') : -7.0, ('F', 'T') : -2.0, ('D', 'D') : 12.0,
('Q', 'T') : -3.0, ('W', 'G') : -5.0, ('Q', 'P') : -3.0, ('W', 'C') : -4.0,
('W', 'K') : -8.0, ('H', 'Q') :  2.0, ('Q', 'D') :  0.0, ('W', 'W') : 11.0,
('V', 'L') :  1.0, ('L', 'G') : -2.0, ('W', 'S') : -5.0, ('L', 'K') : -7.0,
('N', 'P') : -4.0, ('H', 'E') : -1.0, ('N', 'T') : -1.0, ('H', 'A') : -3.0,
('Y', 'L') : -2.0, ('Y', 'H') :  3.0, ('G', 'S') :  1.0, ('Y', 'D') : -4.0,
('V', 'Q') : -3.0, ('L', 'T') : -1.0, ('G', 'G') :  9.0, ('G', 'C') : -2.0,
('E', 'N') :  0.0, ('Y', 'T') : -3.0, ('Y', 'P') : -5.0, ('R', 'N') : -3.0,
('V', 'D') : -5.0, ('K', 'R') : -1.0, ('V', 'H') : -5.0, ('I', 'Q') : -3.0,
('V', 'P') : -4.0, ('M', 'C') : -2.0, ('K', 'N') : -2.0, ('V', 'T') :  0.0,
('M', 'G') : -1.0, ('T', 'S') :  1.0, ('I', 'E') : -5.0, ('M', 'K') : -6.0,
('I', 'A') :  0.0, ('N', 'N') : 11.0, ('R', 'R') :  9.0, ('F', 'P') : -5.0,
('L', 'I') :  2.0, ('I', 'S') : -2.0, ('D', 'S') : -4.0, ('L', 'S') : -2.0,
('I', 'H') : -5.0, ('F', 'D') : -5.0, ('D', 'C') : -7.0, ('F', 'H') : -2.0,
('D', 'G') : -2.0, ('F', 'L') :  1.0
}

chars = [aa[0] for aa in matrix] + [aa[1] for aa in matrix]
chars = sorted(list(set(chars)))

output = " ,%s\n" % ",".join(chars)

for i, aa1 in enumerate(chars):
    output += aa1
    for j, aa2 in enumerate(chars):
        if (aa1, aa2) in matrix:
            output += "," + str(matrix[(aa1, aa2)])
        else:
            output += "," + str(matrix[(aa2, aa1)])
    output += "\n"
print(output)

```
> 2 votes

---
Tags: phylogenetics, phylogeny, matrix, substitution, amino-acids
---


# Question
Title: How to merge multiple list into a single dataframe replacing common values with 1 -python
I have multiple lists on python:

```
a = ["house","garden", "living room", "dog","cat"]

b= ["cat","dog", "chicken"]

c=["house", "garden","bathroom"]

```

I'd like to create the following dataframe:

| a | b | c |
| --- | --- | --- |
| house | 0 | 1 |
| garden | 0 | 1 |
| living room | 0 | 0 |
| dog | 1 | 0 |
| cat | 1 | 0 |

Thanks for the help!

# Answer
Here you go, just add column 'C' and you are sorted. Just to note, you must be seeking a `pandas` solution because `merge` is a `pandas` command and you're asking for a `Python` solution.

```
import pandas as pd

a = ["house","garden", "living room", "dog","cat"]
b= ["cat","dog", "chicken"]
df = pd.DataFrame(a, columns = ['a'])
df2 = pd.DataFrame(b, columns = ['b'])
dfa = df['a'].value_counts()
dfa.columns = ['a']
dfb = df2['b'].value_counts()
dfb.columns = ['b']
dfa = dfa.to_frame()
dfb = dfb.to_frame()
df3 = dfa.join(dfb).replace(np.nan, 0).astype(int)
print (df3)

```

Output

```
             a  b
house        1  0
garden       1  0
living room  1  0
dog          1  1
cat          1  1

```

and if you want to remove the 'a' column

```
print (df3.drop('a', axis=1))

             b
house        0
garden       0
living room  0
dog          1
cat          1

```

---

**Notes** The key command here is `value_counts` and makes a frequency plot. Its output is a `Series` rather than a `DataFrame` so it needs converting.

If you want to keep all inputs the commands is

```
df3 = pd.concat([dfa, dfb]).replace(np.nan, 0).astype(int)

```

---

Alternatively,

```
a = ["house","garden", "living room", "dog","cat"]
b= ["cat","dog", "chicken"]
df = pd.DataFrame(a).value_counts().to_frame('a')
df2 = pd.DataFrame(b).value_counts().to_frame('b')
df3 = df.join(df2).replace(np.nan, 0).astype(int).rename_axis(None)
print (df3.drop('a', axis=1))

```
> 1 votes

# Answer
I'm aware that this asks specifically for a python solution, but I'd like to add an R answer as well, just in case someone using R has a similar problem:

```
> a <- c("house","garden", "living room", "dog","cat")
> b <- c("cat","dog", "chicken")
> c <- c("house", "garden","bathroom")

> (result <- data.frame(row.names=a, a=a %in% a, b=a %in% b, c=a %in% c))
               a     b     c
house       TRUE FALSE  TRUE
garden      TRUE FALSE  TRUE
living room TRUE FALSE FALSE
dog         TRUE  TRUE FALSE
cat         TRUE  TRUE FALSE

```

or for strictly what you asked for, the TRUE/FALSE values can be encouraged into numbers by adding zero:

```
> (result <- data.frame(a=a, b=(a %in% b) + 0, c=(a %in% c) + 0))
            a b c
1       house 0 1
2      garden 0 1
3 living room 0 0
4         dog 1 0
5         cat 1 0

```
> 1 votes

# Answer
Both existing answers seem to work to my eye. An alternative approach:

```
import pandas as pd

a = ["house","garden", "living room", "dog","cat"]
b = ["cat","dog", "chicken"]
c = ["house", "garden","bathroom", "unique"]

```

First we generate all options available. This way your 'a' list is not restricting the options. We use sets to get unique values.

```
options = list(set(a + b + c))
wanted = pd.DataFrame({'options' : options} )

```

Next we loop over all the lists. We use tuples to assign both name and variable.

```
for name, _ in [('a',a),('b',b),('c',c)]:
    wanted.loc[wanted['options'].isin(_), name ] = 1

```

Finally we fill the missing values with zeroes.

```
wanted = wanted.fillna('0')

```

In case you have large number of lists, you can iterate over ascii\_lowercase letters:

```
my_tuples = ()
from string import ascii_lowercase

for _ in ascii_lowercase[:3]: #<-magic number for number of lists here
     my_tuples = my_tuples + (str(_), eval(_))

```

Then change for loop as:

```
for name, _ in my_tuples:
     ....

```
> 1 votes

---
Tags: python, phylogenetics, merge, pandas, dataframe
---


# Question
Title: How to split a Seurat cluster in several subclusters?
I've analyzed my scRNA-seq data and have a couple of Seurat clusters that show more than one cell type in each cluster. (for example, cluster 9 shows both NK and CD4 cells) How can I split a cluster into subclusters, based on the gene expression?

# Answer
> 1 votes
Increase the clustering `resolution` parameter to generate more (smaller) clusters, see FindClusters in the Seurat docs. Whether or not this will neatly, split your clusters into subclusters depends on your data, but normally one can easily separate CD4 and NK cells from PBMCs. See also the Clustree approach for determining the optimal resolution.

Alternatively/additionally use a cell type annotation approach that does not require clustering, e.g. Seurat V4 permits reference mapping and there are many other methods too

# Answer
> 1 votes
I've done sub-clustering a few times on my Seurat data sets. The approach I take is to subset the clusters that need to be clustered (i.e. using `subset`), carry out a clustering of only those cells, then transfer the subcluster labels back to the original dataset. Here's some rough code, which will need to be modified for your specific situation and code preferences:

```
data.sc[["main.cluster"]] <- Idents(data.sc);
clusterToSubset <- "B cells";
data.sub <- subset(data.sc, subset = main.cluster == clusterToSubset);
## ** [Cluster data.sub] **
Idents(data.sc) <- "main.cluster";
## This should work as long as cell order (i.e. column names) isn't changed
Idents(data.sc, cells=colnames(data.sc) %in% colnames(data.sub)) <-
    paste0("Sub_", Idents(data.sub));
data.sc[["merged.cluster"]] <- Idents(data.sc);

```

This approach will work with multiple sub-clusters as well by changing the label that is pasted onto the sub-cluster idents:

```
data.sc[["main.cluster"]] <- Idents(data.sc);
clusterToSubset1 <- "B cells";
data.sub.1 <- subset(data.sc, subset = main.cluster == clusterToSubset1);
## ** [Cluster data.sub.1] **
clusterToSubset2 <- "T cells";
data.sub.2 <- subset(data.sc, subset = main.cluster == clusterToSubset2);
## ** [Cluster data.sub.2] **
Idents(data.sc) <- "main.cluster";
## This should work as long as cell order (i.e. column names) isn't changed
Idents(data.sc, cells=colnames(data.sc) %in% colnames(data.sub.1)) <-
    paste0("Sub1_", Idents(data.sub.1));
Idents(data.sc, cells=colnames(data.sc) %in% colnames(data.sub.2)) <-
    paste0("Sub2_", Idents(data.sub.2));
data.sc[["merged.cluster"]] <- Idents(data.sc);

```

# Answer
> 0 votes
`Idents(data.sc, cells=colnames(data.sc) %in% colnames(data.sub.1))` gives the same output as `Idents(data.sc, cells=colnames(data.sc) %in% colnames(data.sub.2))`.

Both are returning the cellIDs of all cells.

The correct code to generate a list of colnames of data.sc is as follows:

```
Idents(data.sc, cells=colnames(data.sc)[colnames(data.sc) %in% colnames(data.sub.1)]) <- paste0("sub1_", Idents(data.sub.1))

```

However, you can directly add the new cluster names to a metadata column in the Seurat object using the following code:

## Set up new metadata column to store subcluster IDs later ##
data.sc$subcluster <- Idents(data.sc)

## subclustering steps (from a previous comment)##
data.sub.1 <- subset(data.sc, subset = seurat_clusters == 1)
data.sub.1 <- FindVariableFeatures(data.sub.1)
data.sub.1 <- ScaleData(data.sub.1)
data.sub.1 <- RunPCA(data.sub.1, verbose = FALSE)
data.sub.1 <- RunUMAP(data.sub.1, reduction = "pca", dims = 1:30)
data.sub.1 <- FindNeighbors(data.sub.1, reduction = "pca", dims = 1:30)
data.sub.1 <- FindClusters(data.sub.1, resolution = 0.5)

## merging the subclusster names back to original clusters ##
data.sc$subcluster[colnames(dta.sc)[colnames(data.sc) %in% colnames(data.sub.1)]] <- paste0("sub1_", data.sub.1$seurat_clusters)
```

---
Tags: scrnaseq, seurat, single-cell, gene-expression, clustering
---


# Question
Title: How can I convert my tree to the correct nexus format?
I would like to use the program BayesTraitV4. I have my phylogenetic tree in NEXUS format but the program does not accept it. I checked the example data and while they are also in NEXUS format, they seemed different. I will provide a part of my data and the example data respectively.

My tree:

```
#NEXUS
begin taxa;
    dimensions ntax=85;
    taxlabels
    LevantinaFC4779
    Assyriella_naegelei
    Ceu_LAG_2
    CeuMist2
    CeuMist1
    Ceu_PAN_2
    Ceu_PAN_1
    Ceu_PAN_3
    Ceu_PAN_4
    CcMath1
    CcNav1
    CcNav2
    Cc_Mal_1
    Cc_KalR2_2
    CcRod1
    CcRod2
    CcRod3
    CcRod4
    CcHr1
    CeST6_1
    CeST6_4
    CeST6_5
    CeAA7_1
    CeST6_3
    CeST6_2
    Cg_MelR7_1
    Cg_MelR3_1
    Cg_MelR6_2
    Cg_MelR7_2
    Cg_MelR6_1
    Cg_MelR7_3
    Cg_MelR5_1
    Cg_MelR4_2
    Cg_MelR1_1
    Cg_MelR2_1
    Cg_MelR4_1
    Cg_MelR7_4
    Cg_MelR5_2
    ChTM4_2
    ChTM5_1
    ChTM4_3
    ChTM4_1
    Ch_MEN_5
    Ch_MEN_13
    Ch_MEN_8
    Ch_LEV_3
    Ch_LEV_5
    Ch_LEV_4
    Ch_MEN_6
    Ch_MEN_10
    Ch_LEV_2
    Ch_MEN_4
    Ch_MEN_14
    Ch_LEV_1
    Ch_MEN_11
    Ch_MEN_3
    Ch_MEN_9
    Ci_MSP_3
    Ci_MSP_2
    Ci_MSP_1
    Ci_DIA_1
    Ci_DIA_2
    CiPLAN1
    CiPLAN3
    CiPLAN2
    Ci_MEN_1
    Ci_MEN_2
    Ci_MEN_12
    CpGravR6_1
    CpGravR3
    CpGravR4
    CpGio_2
    CpGravR5_2
    CpGravR2
    CpGravR7
    CpGio_1
    CpGio_3
    CpGravR1
    CpGravR6_2
    CpGravR5_1
    CpThe_1
    CpSka_1
    Cn_IG_2
    Cn_IG_3
    Cn_IG_1
;
end;

begin trees;
    tree TREE1 = [&R] ((LevantinaFC4779[&rate_range={7.888601174550756E-4,0.2727277039812744},height_95%_HPD={0.0,1.5711876244495215E-12},length_range={0.6843282088933,6.197383612374},height_median=6.004086117172847E-13,length_95%_HPD={1.080900382604,3.875114824373},height=6.614724410322848E-13,rate=0.049900260878735954,height_range={0.0,5.799805080641818E-12},rate_median=0.045084434131968686,length=2.411967820837161,length_median=2.3472043999879997,rate_95%_HPD={0.0108358382756707,0.10141567250363237}]:2.4119678208371664,Assyriella_naegelei[&rate_range={2.4004562386929353E-4,0.19480057878691465},height_95%_HPD={0.0,1.5711876244495215E-12},length_range={0.6843282088933,6.197383612374},height_median=6.004086117172847E-13,length_95%_HPD={1.080900382604,3.875114824373},height=6.614724410322848E-13,rate=0.03027613080757391,height_range={0.0,5.799805080641818E-

```

Example tree:

```
#NEXUS
Begin trees;
    translate
        1   Opossum,
        2   Diprotodontian,
        3   Sloth,
        4   Armadillo,
        5   Anteater,
        6   Hedgehog,
        7   Mole,
        8   Shrew,
        9   Tenrecid,
        10  Golden_Mole,
        11  Sirenian,
        12  Hyrax,
        13  Elephant,
        15  Lo_Ear_Ele_shrew,
        17  Sciurid,
        18  Mouse,
        19  Rat,
        20  Hystricid,
        21  Caviomorph,
        22  Rabbit,
        23  Pika,
        24  Flying_Lemur,
        25  Tree_Shrew,
        26  Strepsirrhine,
        28  Phyllostomid,
        29  Free_tailed_bat,
        30  False_vampire_bat,
        31  Flying_Fox,
        32  Rousette_Fruitbat,
        33  Whale,
        35  Hippo,
        36  Llama,
        37  Ruminant,
        38  Pig,
        39  Horse,
        40  Rhino,
        41  Tapir,
        42  Cat,
        43  Caniform,
        44  Pangolin;
        tree tree.164000.208997.409491 = ((1:0.07628,2:0.065151):0.1546885,(((15:0.133631,(9:0.157087,10:0.077841):0.007403):0.005959,(11:0.0406,(12:0.092769,13:0.052284):0.001099):0.01453):0.029813,((4:0.05807,(3:0.054688,5:0.062757):0.008238):0.029441,(((26:0.083487,(24:0.064103,25:0.110063):0.003501):0.000775,((22:0.058663,23:0.116805):0.045375,(17:0.0898,((18:0.039267,19:0.043896):0.125808,(20:0.083297,21:0.089139):0.067194):0.006543):0.011062):0.005229):0.00803,((7:0.090648,(6:0.150383,8:0.133021):0.006503):0.01198,((36:0.062932,(38:0.066525,(37:0.076473,(35:0.039222,33:0.031429):0.00357):0.008228):0.005968):0.026714,((44:0.081521,(42:0.046822,43:0.059428):0.024607):0.004257,((39:0.046001,(40:0.029049,41:0.029533):0.003965):0.022397,((28:0.093611,29:0.050095):0.013432,(30:0.061026,(31:0.018601,32:0.018701):0.036144):0.004989):0.015808):0.001216):0.000872):0.002868):0.009567):0.010713):0.00234):0.1546885);
        tree tree.165000.208997.199794 = ((1:0.07628,2:0.065151):0.1546885,(((15:0.133631,(9:0.157087,10:0.077841):0.007403):0.005959,(11:0.0406,(12:0.092769,13:0.052284):0.000485):0.01453):0.029813,((4:0.05807,(3:0.054688,5:0.062757):0.008238):0.029441,(((26:0.083487,(24:0.064103,25:0.110063):0.002826):0.000775,((22:0.058663,23:0.116805):0.045375,(17:0.0898,((18:0.039267,19:0.043896):0.125808,(20:0.083297,21:0.089139):0.067194):0.006543):0.011062):0.003736):0.00803,((7:0.090648,(6:0.150383,8:0.133021):0.010132):0.01198,((44:0.084486,(42:0.046822,43:0.059428):0.024607):0.004257,((36:0.062932,(38:0.066525,(37:0.076473,(35:0.039222,33:0.031429):0.004674):0.008228):0.005968):0.026714,((39:0.046001,(40:0.029049,41:0.029533):0.003965):0.022397,((28:0.093611,29:0.050095):0.008576,(30:0.061026,(31:0.018601,32:0.018701):0.036144):0.004989):0.015808):0.001216):0.000872):0.002868):0.009567):0.010713):0.00234):0.1546885);
;

end;

```

Why are these formats different and how can I make my file to look like the example one? Keep in mind the files are not presented in their entirety because I'm limited with the amount of characters.

Edit: I will be providing my tree in this link https://drive.google.com/file/d/1lYmzW9htwaclC\_yYs3fQom-Z3zbjk8wH/view?usp=sharing

# Answer
> 2 votes
Okay ... this will probably work

```
#NEXUS
begin trees;
    tree TREE1 = [&R] ((LevantinaFC4779:2.411968,Assyriella_naegelei:2.411968):5.752052,((Cn_IG_1:0.126901,(Cn_IG_2:0.075971,Cn_IG_3:0.075971):0.05093):7.4805,(((CpThe_1:0.122189,CpSka_1:0.122189):2.663216,((CpGio_2:0.137489,(CpGravR6_1:0.087156,(CpGravR3:0.038234,CpGravR4:0.038234):0.048922):0.050334):0.257308,(CpGravR5_1:0.28159,((CpGravR7:0.099712,(CpGravR5_2:0.020255,CpGravR2:0.020255):0.079458):0.15617,(CpGravR6_2:0.109222,(CpGio_1:0.068711,(CpGio_3:0.049345,CpGravR1:0.049345):0.019367):0.040511):0.14666):0.025708):0.113207):2.390607):1.772627,(((Ci_MEN_12:0.072368,(Ci_MEN_1:0.022725,Ci_MEN_2:0.022725):0.049643):1.78894,((CiPLAN2:0.253487,(CiPLAN1:0.067509,CiPLAN3:0.067509):0.185978):0.386279,((Ci_MSP_3:0.061614,Ci_MSP_2:0.061614):0.231268,(Ci_MSP_1:0.113777,(Ci_DIA_1:0.047073,Ci_DIA_2:0.047073):0.066704):0.179105):0.346885):1.221542):1.580745,(((Ceu_LAG_2:1.580984,(CeuMist2:0.198982,CeuMist1:0.198982):1.382002):0.341226,((Ceu_PAN_2:0.046898,Ceu_PAN_1:0.046898):0.20368,(Ceu_PAN_3:0.018814,Ceu_PAN_4:0.018814):0.231764):1.671632):1.226475,((CcHr1:1.247104,((CcRod1:0.403751,(CcRod4:0.276104,(CcRod2:0.158977,CcRod3:0.158977):0.117127):0.127647):0.308704,((Cc_Mal_1:0.138232,Cc_KalR2_2:0.138232):0.395435,(CcMath1:0.338389,(CcNav1:0.091762,CcNav2:0.091762):0.246627):0.195278):0.178788):0.534649):1.609976,(((CeST6_5:0.056918,(CeST6_1:0.030107,CeST6_4:0.030107):0.026812):0.733078,((CeST6_2:0.143014,(CeAA7_1:0.097814,CeST6_3:0.097814):0.0452):0.349431,((Cg_MelR7_1:0.136715,(Cg_MelR3_1:0.088065,(Cg_MelR6_2:0.047987,(Cg_MelR7_2:0.019268,Cg_MelR6_1:0.019268):0.028719):0.040078):0.04865):0.074941,(Cg_MelR5_2:0.143401,(Cg_MelR7_3:0.119131,(Cg_MelR5_1:0.090065,((Cg_MelR4_1:0.020938,Cg_MelR7_4:0.020938):0.051613,(Cg_MelR4_2:0.052913,(Cg_MelR1_1:0.021387,Cg_MelR2_1:0.021387):0.031526):0.019639):0.017514):0.029066):0.02427):0.068255):0.280789):0.297551):0.793423,(((ChTM4_2:0.049871,ChTM5_1:0.049871):0.111934,(ChTM4_3:0.023134,ChTM4_1:0.023134):0.138672):0.228306,(Ch_MEN_9:0.174548,((Ch_MEN_11:0.021817,Ch_MEN_3:0.021817):0.10954,((Ch_MEN_8:0.057141,(Ch_MEN_5:0.022973,Ch_MEN_13:0.022973):0.034168):0.053065,((Ch_LEV_3:0.031056,Ch_LEV_5:0.031056):0.061225,((Ch_LEV_4:0.029359,(Ch_MEN_6:0.017456,Ch_MEN_10:0.017456):0.011904):0.034407,((Ch_LEV_2:0.017911,Ch_MEN_4:0.017911):0.040741,(Ch_MEN_14:0.03419,Ch_LEV_1:0.03419):0.024462):0.005114):0.028515):0.017925):0.021151):0.043191):0.215563):1.193307):1.27366):0.291605):0.293368):1.115978):3.04937):0.556619);
end;

```

Alternatively, I'm not sure about the root, just remove `[&R]` in the above.

Finally there could be an issue that you haven't obtained the sample directly from the raw Beast output, instead you've summarised it. It might complain about that. Thus it might need all those trees from the raw Beast output. Thats not a format issue it's a data input issue. BayesTraits might be perform the calculation against the total MCMCMC sample minus the burn-in.

---

from comments

Read it into Rs `ape::read.nexus()` and then write via `phytools::writeNexus()`. This should produce your translate table. Biopython recognises this format but doesn't write to it.

The closest Biopython can achieve is a taxablock here

```
#NEXUS
Begin Taxa;
 Dimensions NTax=85;
 TaxLabels LevantinaFC4779 Assyriella_naegelei Cn_IG_1 Cn_IG_2 Cn_IG_3 CpThe_1 CpSka_1 CpGio_2 CpGravR6_1 CpGravR3 CpGravR4 CpGravR5_1 CpGravR7 CpGravR5_2 CpGravR2 CpGravR6_2 CpGio_1 CpGio_3 CpGravR1 Ci_MEN_12 Ci_MEN_1 Ci_MEN_2 CiPLAN2 CiPLAN1 CiPLAN3 Ci_MSP_3 Ci_MSP_2 Ci_MSP_1 Ci_DIA_1 Ci_DIA_2 Ceu_LAG_2 CeuMist2 CeuMist1 Ceu_PAN_2 Ceu_PAN_1 Ceu_PAN_3 Ceu_PAN_4 CcHr1 CcRod1 CcRod4 CcRod2 CcRod3 Cc_Mal_1 Cc_KalR2_2 CcMath1 CcNav1 CcNav2 CeST6_5 CeST6_1 CeST6_4 CeST6_2 CeAA7_1 CeST6_3 Cg_MelR7_1 Cg_MelR3_1 Cg_MelR6_2 Cg_MelR7_2 Cg_MelR6_1 Cg_MelR5_2 Cg_MelR7_3 Cg_MelR5_1 Cg_MelR4_1 Cg_MelR7_4 Cg_MelR4_2 Cg_MelR1_1 Cg_MelR2_1 ChTM4_2 ChTM5_1 ChTM4_3 ChTM4_1 Ch_MEN_9 Ch_MEN_11 Ch_MEN_3 Ch_MEN_8 Ch_MEN_5 Ch_MEN_13 Ch_LEV_3 Ch_LEV_5 Ch_LEV_4 Ch_MEN_6 Ch_MEN_10 Ch_LEV_2 Ch_MEN_4 Ch_MEN_14 Ch_LEV_1;
End;
Begin Trees;
 Tree tree1=((LevantinaFC4779:2.41197,Assyriella_naegelei:2.41197):5.75205,((Cn_IG_1:0.12690,(Cn_IG_2:0.07597,Cn_IG_3:0.07597):0.05093):7.48050,(((CpThe_1:0.12219,CpSka_1:0.12219):2.66322,((CpGio_2:0.13749,(CpGravR6_1:0.08716,(CpGravR3:0.03823,CpGravR4:0.03823):0.04892):0.05033):0.25731,(CpGravR5_1:0.28159,((CpGravR7:0.09971,(CpGravR5_2:0.02025,CpGravR2:0.02025):0.07946):0.15617,(CpGravR6_2:0.10922,(CpGio_1:0.06871,(CpGio_3:0.04934,CpGravR1:0.04934):0.01937):0.04051):0.14666):0.02571):0.11321):2.39061):1.77263,(((Ci_MEN_12:0.07237,(Ci_MEN_1:0.02272,Ci_MEN_2:0.02272):0.04964):1.78894,((CiPLAN2:0.25349,(CiPLAN1:0.06751,CiPLAN3:0.06751):0.18598):0.38628,((Ci_MSP_3:0.06161,Ci_MSP_2:0.06161):0.23127,(Ci_MSP_1:0.11378,(Ci_DIA_1:0.04707,Ci_DIA_2:0.04707):0.06670):0.17910):0.34688):1.22154):1.58075,(((Ceu_LAG_2:1.58098,(CeuMist2:0.19898,CeuMist1:0.19898):1.38200):0.34123,((Ceu_PAN_2:0.04690,Ceu_PAN_1:0.04690):0.20368,(Ceu_PAN_3:0.01881,Ceu_PAN_4:0.01881):0.23176):1.67163):1.22647,((CcHr1:1.24710,((CcRod1:0.40375,(CcRod4:0.27610,(CcRod2:0.15898,CcRod3:0.15898):0.11713):0.12765):0.30870,((Cc_Mal_1:0.13823,Cc_KalR2_2:0.13823):0.39543,(CcMath1:0.33839,(CcNav1:0.09176,CcNav2:0.09176):0.24663):0.19528):0.17879):0.53465):1.60998,(((CeST6_5:0.05692,(CeST6_1:0.03011,CeST6_4:0.03011):0.02681):0.73308,((CeST6_2:0.14301,(CeAA7_1:0.09781,CeST6_3:0.09781):0.04520):0.34943,((Cg_MelR7_1:0.13672,(Cg_MelR3_1:0.08807,(Cg_MelR6_2:0.04799,(Cg_MelR7_2:0.01927,Cg_MelR6_1:0.01927):0.02872):0.04008):0.04865):0.07494,(Cg_MelR5_2:0.14340,(Cg_MelR7_3:0.11913,(Cg_MelR5_1:0.09007,((Cg_MelR4_1:0.02094,Cg_MelR7_4:0.02094):0.05161,(Cg_MelR4_2:0.05291,(Cg_MelR1_1:0.02139,Cg_MelR2_1:0.02139):0.03153):0.01964):0.01751):0.02907):0.02427):0.06825):0.28079):0.29755):0.79342,(((ChTM4_2:0.04987,ChTM5_1:0.04987):0.11193,(ChTM4_3:0.02313,ChTM4_1:0.02313):0.13867):0.22831,(Ch_MEN_9:0.17455,((Ch_MEN_11:0.02182,Ch_MEN_3:0.02182):0.10954,((Ch_MEN_8:0.05714,(Ch_MEN_5:0.02297,Ch_MEN_13:0.02297):0.03417):0.05307,((Ch_LEV_3:0.03106,Ch_LEV_5:0.03106):0.06123,((Ch_LEV_4:0.02936,(Ch_MEN_6:0.01746,Ch_MEN_10:0.01746):0.01190):0.03441,((Ch_LEV_2:0.01791,Ch_MEN_4:0.01791):0.04074,(Ch_MEN_14:0.03419,Ch_LEV_1:0.03419):0.02446):0.00511):0.02851):0.01792):0.02115):0.04319):0.21556):1.19331):
1.27366):0.29161):0.29337):1.11598):3.04937):0.55662):0.00000;
End;

```

The solution could be coded from the taxa block if the above doesn't work, but easier to use **R** `ape` and `phytools`

---
Tags: phylogenetics, phylogeny, file-formats
---


# Question
Title: How to measure the total size of a fastq file in base pairs?
Or Kbps/Gbps. It feels like it should be conceptually very simple, but I can't seem to figure out the right combination of keywords to find it via my search engine. Help would be appreciated!

I have BBMAP, SRAtoolkit and MEGAHIT already installed, and also use bash. I'd be very happy if this can be answered with software that I already have, but if not that's perfectly fine.

# Answer
> 5 votes
I've been using this:

`cat file.fastq | paste - - - - | cut -f 2 | tr -d '\n' | wc -c`

Explanation :

* `paste - - - -` : print four consecutive lines in one row (tab delimited), to merge the info for each read
* `cut -f2` : print only the second column, to access the sequence after the paste
* `wc -c` : count the characters
* `tr -d '\n'`: to remove from count the eventual newline characters

(a tip for your googling: try search for "counting number of bases in fastq file")

# Answer
> 1 votes
The number of bases in a fastq file can be counted in bash with awk and wc

```
awk 'NR % 4 == 0' ORS="" fastqfile|wc -m

```

The awk code prints every fourth line (which is actually the quality scores, but that doesn't matter here). wc -m returns the total number of characters.

Edited to add ORS="" to prevent count of newlines (the double quotes are not strictly required)

# Answer
> 1 votes
**Using Perl:**
Print the sequence lines (line number 2, 6, 10, etc). Remove the newlines with `chomp`. Count the bytes (here, bases) using `wc -c`:

```
perl -ne 'if ( $. % 4 == 2 ) { chomp; print; }' | wc -c

```

# Answer
> 1 votes
Using a fastq parser in Python (for instance the one provided by `mappy`), it is rather straightforward to compute the sum of the lengths of the sequences present in a fastq file:

```
$ python3 -m pip install mappy  # Only if mappy is not installed yet
$ python3 -c 'from mappy import fastx_read; print(sum(len(seq) for (_, seq, *_) in fastx_read("SRR077487_2.filt.fastq.gz")))'
2386161200

```

`python3 -c '<insert some python code here>'` uses the `-c` option of `python3` that tells Python to execute the provided code (instead of starting an interactive interpreter, or running code present in a file). Here the code has two steps (separated by a `;`, in a standard script, we could use a new line instead): first making the `fastx_read` function available (because it is part of an optional module), second displaying the total length of the sequences.

`fastx_read("path_to_a_file")` generates (name, sequence, quality) triplets when a fastq file is provided, or (name, sequence) pairs when a fasta file is provided.

The `(_, seq, *_)` syntax is a form of "tuple unpacking" where we store the second element (the sequence) in a `seq` variable. The `_` in is just a way to ignore the first element in the triplet/pair. The `*_` is a way to ignore whatever is after the second element in the triplet/pair.

`len(seq) for (_, seq, *_) in fastx_read("SRR077487_2.filt.fastq.gz")` is an expression representing the lengths (obtained with the predefined `len` function) of the sequences present in the file `"SRR077487_2.filt.fastq.gz"`.

So the code uses the predefined `sum` function to compute the sum of values that are the lengths of the sequences extracted from the fastq file, and this is displayed using the `print` function.

This should also work with either fasta or fastq format, and with or without gz compression.

I tested more parsers on a slightly different problem here: https://bioinformatics.stackexchange.com/a/380/292 and you can find inspiration with the other approaches proposed in the other answers there.

# Answer
> 1 votes
`seqkit stats fastq.gz`. Get SeqKit from https://bioinf.shenwei.me/seqkit/usage/

# Answer
> 0 votes
With bbmap:

`reformat.sh </path/to/fastq>`

---
Tags: fastq, bash
---


# Question
Title: Is 16S sequencing cheaper than whole-genome sequencing for taxonomic assignment?
*This is adapted from a question posted on reddit*

In my work, I occasionally get asked about metagenomic sequencing, which (after further questioning) I discover means that the person wants to work out the taxonomic composition of a sample.

In the past, people have gone straight to 16S sequencing for this, and various manufacturers have created targeted 16S kits for this purpose. But does it still make sense to do that? Is 16S still cheaper and better than other options?

# Answer
> 2 votes
When comparing anything targeted to whole-genome sequencing for the purposes of assembly, the targeted approach will always be cheaper, because the total length that needs to be assembled is shorter. But that's like comparing a video to a photograph.

If you're just looking to determine the phylogenetic composition of a bacterial sample, then there's no reason why targeted 16S sequencing should be cheaper than whole-genome shotgun sequencing, and the whole genome approach may end up better due to the highly conserved nature of the 16S sequence, and the abundance of horizontal gene transfer.

Previously, 16S might have been better because the taxonomic databases were more complete, but now that there are so many microbial genomes available, the whole genome databases are comprehensive enough for most purposes.

@terdon linked to this 2008 paper, in which they compare random sequence reads to 16S for estimating the biodiversity of a metagenomic library, finding that 16S wasn't worth it and they could get perfectly good results with random reads.

In a situation where samples are abundant and no amplification is required, PCR can be excluded from whole-genome sequencing, making it the cheaper option at the sample preparation side of things.

Just to put numbers to things to demonstrate how cheap whole-genome sequencing can be, here are the big consumable costs for nanopore (prices in USD):

* RBK114.24 rapid barcoding kit\* - \\$700 (24 samples, 12 runs)
* Flongle flow cell - \\$810 (12 flow cells)

So the marginal cost of one run is \\$68 for the flow cell, and \\$60 for the sequencing kit, processing 24 samples at a cost of about $11 per sample (rounding up).

\[FWIW, ONT also provide a rapid PCR barcoding kit for the same cost (i.e. \\$700, but polymerase is not included), which can be used together with PCR for low-yield samples.\]

Flongle flow cells will usually generate \>200 Mb per run (up to about 2 Gb in rare cases), so that's about 8 Mb of sequence per sample, which should be sufficient to do proportional counting down to about 1% of the sample (assuming \>1000 reads per sample).

* There is a 96-barcode kit available for \\$990, but I don't think that would have enough coverage to be useful when running samples on a Flongle flow cell.

---
Tags: phylogenetics, metagenome, hts
---


# Question
Title: Does phasing improve imputation quality?
If so, by about how much? I want to know the downsides of skipping the phasing step on human microarray and WGS data.

# Answer
> 5 votes
Yes. Phasing data leads to a better imputation accuracy. More specifically:

* Improves allele matching: Correct phased data ensures that alleles are matched correctly between the reference panel and the target dataset.
* Using phased data, the computational time of imputation is significantly reduced.
* Having haplotypes (phased data) helps imputing variants that were not genotyped (missing variants/genotypes).

This paper based on the UK10K data has some actual numbers describing imputation accuracy of phased vs unphased data.

---
Tags: microarray, wgs, imputation, phasing
---


# Question
Title: What does the different sequences represent?
I am using this package nsdpy to download genome sequences from **NCBI nucleotide database**.
Specifically I am interested in the **whole mitochondrial genome** of different species, here I will use a subset of species just as an example.

My Query is:
`"{mitochondrion[Title]" -i -t -L subset_species_list.txt`
Where :

* `-i` Adds the taxonomic information to the information lines of the sequences written in the output files (not relevant for this example)
* `-t` option produces a tsv file with the following fields: Taxon name, SequenceID, TaxID, Taxonomic lineage, Sequence length, Sequence.
* `subset_species_list.txt` is a list of species, below you can see its content

```
Xenopus laevis
Gasterosteus aculeatus
Chironomus riparius
Rasbora heteromorpha
Planktolyngbya limnetica

```

The output file is quite long and includes the actual sequences. Here I will post the same file where I dropped the column for the sequence and the lineage to keep it compact

My aim is to obtain phylogenetic distances between species.

My questions are:

1. What do the sequences for the same species represent?
2. Which sequences are best suited for an alignment and how do I know that?
3. Should I also look into the orientation of the sequences or are they all on the forward strand?
4. Is it better if I filter for Coding sequences (CDS) instead of whole genomes

Also if you noticed any rookie mistake in my process feel free to point that out. It is the first time for me working with DNA sequences and I really want to learn to correct way to do that.

# Answer
> 0 votes
Welcome! I have some quick basic answers for the questions you're asking so hopefully this helps:

1. Sequences for the same species represent variation/diversity/evolutionary changes. Think of humans; how do different people's DNA sequences represent the sequence of a human? This is especially important for phylogenetic studies.
2. For the most part, all sequences are best suited for alignment. The main thing to think about is the amount of data you have (if it's a lot of GB the alignment will take a lot of computing power), the type of sequnece data you have (simple rule of thumb is to not mix short read with long read sequences), and the tool you're using for alignment.
3. Always a good idea to think about orientation but most tools will account for orientation when aligning sequences.
4. This depends on what you're looking for but more often than not using the whole genomes is best as it lets you look at the whole picture instead of just a small portion but that's definitely a lot of computing power. If you're just looking at the mitochondria just focus on that region

Hope this helps!

---
Tags: sequence-alignment, sequence-analysis
---


title: Summary statistics version problem

text:
# Question
Title: Summary statistics version problem
I am learning GWAS study and PRS for predicting disease. I was aligning my sequencing data with GRCH38 reference panel. However, in the PRS analysis step, I found that most of the summary statistics are developed in GRCH37. It was a bit inconvenient for me because I have to convert vcf/snp file into GRCH37 which may end up losing some snps. I was checking this database link. Is there any database of summary statistics in GRCH38?

# Answer
This is a common issue. Instead of lifting over your genomes, it might be computationally and conceptually easier just to lift over the summary statistics to GRCh37. Alternatively (and preferably), this is **done for you** for every GWAS published in the GWAS catalog. Simply click on the "harmonized" folder after clicking on "download" and you can download whatever build aligns with your genotype data.
> 2 votes

# Answer
You might want to check out this tool which contains some software to liftover your coordinates from one build to another.
> 2 votes

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


# Question
Title: How to associate point mutations to quantitative protein levels
My experiment is based on a quantitative continuous variable levels (B12 vitamin) in a sample of patients retrospectively selected based on values that can take (under our criteria):

* normal
* low levels: under the physiological or normal reference values but not enough to be considered deficiency
* deficiency (which is associated with poor nutrition, gastric disease, etc.). These data are excluded.

We are interested in **low levels**. Our hypothesis is **to prove if mutations in 10 genes can be associated with these physiological low levels.** These physiological low levels are quantified under the format of maximum and minimum value of B12. The genes can take 3 possible categories:

* wt: absence of mutation
* het: half of mutation
* hom: full mutation

Hypothesis:

1. I would like, if possible, to try something like a genetic risk score, **if multiple het or hom states** in multiple genes can be additive. I don't know how to approach this, because there can be multiple permutations, and maybe permutation 1 has a stronger association with low levels than permutation 2.
2. Any ideas on how to carry this out? I don't really know if mutations are additive or opposite. Not sure if some clustering approach can be helpful here, like PCA or euclidean distance, linear discriminant analysis... Any advice is welcome.

My idea is to execute a summary of levels of B12 for each gene and see in a boxplot possible differences, perhaps look for some association (or correlation), but I doubt it.

Any idea or previous input or experiment like this would be tremendously welcome.

# Answer
> 1 votes
Yes, you can run an additive model (i.e., B\*AA = 2\*B\*Aa, aa = B\*AA - 4\*Aa) or a dominance model (B\*AA = B\*Aa). You can test both of these and see which is a better fit to your data.

However, it is practically known already that the model will be at least predominantly additive.

---
Tags: point-mutation, polygenic-risk-score
---


# Question
Title: Approach for alignment-free phylogenetic tree generation using kmers
I have a general question about building phylogenetic trees without alignment:

I want to build a tree using an alignment-free approach by counting kmers. I have found a lot of information in the literature about this type of practice but I'm having trouble finding approaches that use Illumina short read sequences (which is what I have).

In my first attempt, I was able to generate a dendrogram by counting kmers and clustering them, but it did not group the sequences I had in the “correct/expected” way (these are known strains that are part of specific subclasses related to S. cerevisiae - in theory, this kmer clustering approach should’ve been able to group the strains based on their subclass pretty easily).

There are 21 samples in total with 5 specific subclasses. Does anyone know a way to either cluster or simply calculate the distance between kmers to generate a tree of 20+ Illumina short-read samples?

**EDIT:** The main reason for doing this in this way is to eliminate any bias from the reference sequence and maintain all the alleles (most of the samples are tetraploid) so I really want to avoid collapsing the data with any type of preprocessing alignment or assembly.

With that being said I would want to use the whole-genome short-read sequences of these 21 samples to count the kmers across the genome for phylogenetic analysis. If it gets too complicated, I would be okay with assembling the reads first before kmer counting.

I used Kitsune to determine the optimal kmer size for the samples which comes out to 15mer. I was able to count 15mers with all the samples using three methods: pure raw reads, assembling the reads with SPAdes, and downsampling the raw reads to 5x coverage. Using Jellyfish and custom scripts I wrote, which didn't end up being too taxing on my computer/server.

Currently, I have 21 '.tsv' files that list the kmer and the kmer frequency for each individual sample and I have one giant matrix that combines that data together so I have the kmer and the kmer count for each sample in one file:

|  | sample\_1 | sample\_2 | sample\_3 | ... | sample\_21 |
| --- | --- | --- | --- | --- | --- |
| AAAAA | 123432 | 345432 | 56765 | ... | 55677 |
| AAAAG | 1023 | 0 | 350655 | ... | 1 |
| AAAAT | 462 | 741 | 2136 | ... | 121 |
| AAAAC | 342 | 123 | 3453 | ... | 12345 |
| ------------ | ---------- | ---------- | --------------- |  |  |

I have looked briefly at kWIP. Ideally, I would use their distance calculation function to generate a .dist file for my giant matrix - I'll look into what inputs can be used for that function.

# Answer
> 1 votes
A very simple solution is to use a correlation *r* or *r*2. This establishes a pairwise matrix which is then resolved under neighbour-joining (don't use UPGMA or similar). MegaX will import a pairwise matrix and produce a tree in neighbour-joining. You could even produce a bootstrapped solution (that would require custom code) - bit tricky.

**I definitely recommend neighbour-joining for the tree**

MegaX is a good implementation of neighbour-joining (it was originally Nei's idea and that was his package). The reason for this are complex and happy to address this in a separate question.

**Algorithm** Lets say *r* or *r*2 with Sample 1 versus Sample 2 ... the correlation value goes in one half of a triangular matrix with zeros running down the centre of the matrix. You then proceed for every pairwise comparison until you get a triangle with Sample 1, Sample 2, Sample 3, etc ... along the top and Sample 1, Sample 2, Sample 3, etc ... along the side. That will import to MegaX (its a long while since I've used it BTW).

**Distance Methods** If I was performing this I wouldn't use *r* or *r*2: I would use either Jaccard, Euclidean or Cosine (very trendy right now) distances. Cosine is trendy because taps into tensors (Tensorflow and all that). **Scipy** will almost certainly do the distance methods listed here for your data and output a nice triangular matrix.

---

Personally if I was coding this I would absolutely perform a bootstrap, but this is a lot of coding.

**Issues** The issue with distance metrics without bootstrapping is that they give different solutions and there's no way identify which tree is the best tree. Bootstrapping helps a lot here and values \> 80% are considered robust topologies. I've never personally known two distance metrics produce an incongruence for the same data under bootstrapping. The distance matrix issue was why everyone moved (and rigidly stuck with) formal models via maximum likelihood or Bayes - but these require alignments.

**Summary** Think you can go alignment free? Well yeah you can, but remember there's no such thing as a free lunch ;-)

---
Tags: phylogenetics, phylogeny, clustering, k-mer, distance
---


# Question
Title: How to calculate average from a column when consecutive cells are similar in different columns?
I wonder how can I calculate the average of column 7 (WPS\_win) when consecutive cells in column 2, column 3, and column 1 are the same. The file is in bed format (tab-separated format). The original file is very large. The first few lines of the file look like this,

```
chrom   start_pr    end_pr  width_pr    start_win   end_win WPS_win .overlap
chr1    981918  984117  2200    981900  982000  -5  82
chr1    981918  984117  2200    982000  982100  -1  100
chr1    981918  984117  2200    982200  982300  -2  100
chr1    981918  984117  2200    982300  982400  2   100
chr1    981918  984117  2200    982400  982500  -3  100
chr1    981918  984117  2200    982500  982600  -1  100
chr1    981918  984117  2200    982600  982700  -3  100
chr1    981918  984117  2200    982700  982800  1   100
chr1    981918  984117  2200    982800  982900  -5  100
chr1    981918  984117  2200    982900  983000  0   100
chr1    981918  984117  2200    983000  983100  0   100
chr1    981918  984117  2200    983100  983200  -5  100
chr1    981918  984117  2200    983200  983300  -1  100
chr1    981918  984117  2200    983300  983400  -3  100
chr1    981918  984117  2200    983400  983500  1   100
chr1    981918  984117  2200    983500  983600  -2  100
chr1    981918  984117  2200    983600  983700  1   100
chr1    981918  984117  2200    983700  983800  -3  100
chr1    981918  984117  2200    983800  983900  1   100
chr1    981918  984117  2200    983900  984000  -4  100
chr1    981918  984117  2200    984000  984100  -1  100
chr1    981918  984117  2200    984100  984200  -2  17
chr1    999138  1001337 2200    999100  999200  1   62
chr1    999138  1001337 2200    999200  999300  -1  100
chr1    999138  1001337 2200    999700  999800  -1  100
chr1    999138  1001337 2200    999800  999900  -2  100
chr1    999138  1001337 2200    999900  1000000 -2  100
chr1    999138  1001337 2200    1000000 1000100 1   100
chr1    999138  1001337 2200    1000100 1000200 -2  100
chr1    999138  1001337 2200    1000200 1000300 0   100
chr1    999138  1001337 2200    1000300 1000400 -1  100
chr1    999138  1001337 2200    1000400 1000500 -2  100
chr1    999138  1001337 2200    1000500 1000600 -2  100
chr1    999138  1001337 2200    1000600 1000700 2   100
chr1    999138  1001337 2200    1000700 1000800 -1  100
chr1    999138  1001337 2200    1000800 1000900 -2  100
chr1    999973  1002172 2200    999900  1000000 -2  27
chr1    999973  1002172 2200    1000000 1000100 1   100
chr1    999973  1002172 2200    1000100 1000200 -2  100
chr1    999973  1002172 2200    1000200 1000300 0   100
chr1    999973  1002172 2200    1000300 1000400 -1  100
chr1    999973  1002172 2200    1000400 1000500 -2  100
chr1    999973  1002172 2200    1000500 1000600 -2  100
chr1    999973  1002172 2200    1000600 1000700 2   100
chr1    999973  1002172 2200    1000700 1000800 -1  100
chr1    999973  1002172 2200    1000800 1000900 -2  100
chr1    999973  1002172 2200    1001600 1001700 -1  100
chr1    999973  1002172 2200    1001700 1001800 -2  100
chr1    999973  1002172 2200    1001800 1001900 1   100
chr1    999973  1002172 2200    1001900 1002000 -1  100

```

In the above table, the first 22 rows of cols 1, 2, and 3 are the same. Therefore, The expected output would look like this,

```
chrom   start_pr    end_pr  Mean_WPS
chr1    981918  984117  -1.590909091
chr1    999138  1001337 -0.857142857
chr1    999973  1002172 -0.857142857

```

I'm open to both R and Bash. Any suggestion to overcome this issue would be heartily appreciated.

# Answer
```
library(dplyr)
cols_of_interest <- c("chrom", "start_pr", "end_pr")

test_data |>
  group_by(across(all_of(cols_of_interest))) |>
  summarize(Mean_WPS = mean(WPS_win))

```

gives the following with your test data above:

```
# A tibble: 3 × 4
# Groups:   chrom, start_pr [3]
  chrom start_pr  end_pr Mean_WPS
  <chr>    <int>   <int>    <dbl>
1 chr1    981918  984117   -1.59
2 chr1    999138 1001337   -0.857
3 chr1    999973 1002172   -0.857

```

Ref for the code used above
> 2 votes

# Answer
I use `bedtools merge` with a little trick to combine the first three columns into a single one to "mock" each of the unique intervals into one "chromosome" so bedtools easily knows that these are the same:

```
cat example.bed \
| awk 'NR > 1 {OFS="\t"; print $1"--"$2"--"$3"--", 1, 1, $7}' \
| bedtools merge -i - -c 4 -o mean \
| awk 'BEGIN {OFS="\t"; print "chrom", "start_pr", "end_pr", "Mean_WPS"} {split($1,x,"--"); print x[1], x[2], x[3], $4}'

chrom   start_pr        end_pr  Mean_WPS
chr1    981918  984117  -1.590909091
chr1    999138  1001337 -0.8571428571
chr1    999973  1002172 -0.8571428571

```

First `awk` command: Skips first line, combines first three columns into one and prints this together with a single position for start and end plus the score.

`bedtools` command: Merges the identical intervals, taking the mean of the fourth column.

Second `awk` command: Splits the first column again into the former chrom-start-end and puts the average score in column 4, and also adds a header line.

Memory footprint of this should be minimal.
> 0 votes

---
Tags: r, genomics, bash, tsv
---


# Question
Title: detecting DNA of possible infection on cell line
we've done Whole Genome Sequencing with Illumina on DNA extracted from cell lines, where some are suspected to have been infected by something (viral, bacterial, not sure).

I've tried running metaphlan4 (version from github) and got zero hits in all of them.

What other tools could I be trying? Thx.

# Answer
> 3 votes
I've had some success looking for contamination in Illumina and Nanopore reads using kraken2.

Kraken2 has a phylogenetic-based approach, which in the first instance tries to find the most specific taxa that the target read likely comes from (going further up in the tree if there's ambiguity), and also produces a summary table indicating the approximate proportion of reads that come from particular taxa. It can be paired with bracken for additional reliability.

# Answer
> 1 votes
Maybe less advanced but worth mentioning, I've had good results with fastq\_screen. It's simply trying to align a subset of your sequences against various genomes (comes with a few genomes, you can add more).

# Answer
> 1 votes
I didn't manage to try the ones reported by @gringer and @alexlok here, but I did manage to perform an analysis using the OneCodex webapp: https://app.onecodex.com/ and based on the experimental information for the samples being analysed, the results were concordant with the expectation. You can do quite a lot on the "free tier", so at the moment I didn't upgrade to the paid version.

---
Tags: metagenome, whole-genome, dna
---


# Question
Title: Error importing Maximum Credibility Tree from BEAST into R using treeio's read.beast()
I've encountered an issue that I'm hoping someone can help with. I constructed a phylogenetic tree using **BEAST1** and then extracted the Maximum Credibility Tree using treeannotator. I'm trying to import it into **R** for editing, but I run into an error when using the `read.beast()` function from the **treeio** package:

```
MCC_tree <- read.beast("/home/MyPC/Downloads/MCC_beast.tre")
Error in FUN(X[[i]], ...) :
  numbers of left and right parentheses in Newick string not equal

```

Interestingly, when I attempt to visualize the tree using figtree, it works fine. Likewise, importing into R using `read.nexus()` works without issues. However, I want to use `read.beast()` specifically because I need to annotate the "posterior", which I can't do when I import it using `read.nexus()`.

Has anyone encountered a similar issue or have any suggestions on how to resolve this?

Thank you in advance!

# Answer
1. if you read it into FigTree then in the export in Nexus and perform a "save all trees" it should then read into `read.beast()`.
2. Use the beast file output before you performed tree annotator operation.

1 or 2. should be read directly into `read.beast()`. I would try option 2 first.
> 0 votes

---
Tags: r, phylogenetics, phylogeny, bioconductor, ggtree
---


# Question
Title: A function to produce M fragments of DNA with length of I
I want a script in R, which creates `M` fragments of DNA with a length of `I` but with **GC content** of \>70%

I know this gives something but not exactly what I am looking for

```
require(Biostrings)

set = NULL
for (i in 1:4)
  set = c(set, paste(sample(DNA_ALPHABET[1:4], 3, replace=T), collapse=""))
rep(set,2)

```

Any idea please?

# Answer
There are different ways to approach the GC bias, the easiest is probably to decide beforehand how many "A/T" and how many "G/C" we need in each sequence, then choose them randomly.

```
# set global parameters
M <- 7
I <- 15

# for each sequence, choose a proportion of GC
prop_gc <- runif(M,
                 min = .7,
                 max = 1)

# randomly choose bases, ensure we have a number that gives us at least `prop_gc`
bases <- lapply(prop_gc,
                \(.prop_gc){
                  c(
                    sample(c("A","T"),
                           size = floor(I*(1 - .prop_gc)), replace = TRUE),
                    sample(c("G","C"),
                           size = ceiling(I*.prop_gc), replace = TRUE)
                  )
                })

# now we can simply take our bases and combine them (in a random order)
sequences <- sapply(bases,
                    \(.nt) sample(.nt) |> paste(collapse = ""))

sequences
#> [1] "CGGCTTCCGCGCGGGT" "CTACCCACGGCGCCGG" "GGGCGGTGTATGTGCC" "GGCCCGGGTGCCCAGC"
#> [5] "GCGGGCGCGGGCGGGT" "CACCCCCTCCCGGGAC" "GTCGAAGCGCGCCAGC"

```

Created on 2023-10-12 with reprex v2.0.2

Or same thing with a `for` loop (since you gave it as a tag):

```
M <- 7
I <- 15

sequences <- character(length = M)
for(i in seq_len(M)){
  prop_gc <- runif(1,
                   min = .7,
                   max = 1)

  bases <- c(
    sample(c("A","T"),
           size = floor(I * (1 - prop_gc)), replace = TRUE),
    sample(c("G","C"),
           size = ceiling(I * prop_gc), replace = TRUE)
  )

  sequences[[i]] <- sample(bases) |> paste(collapse = "")

}
sequences
#> [1] "CCCCGGTTGCGTTTCC" "CCTCGCTGCGCGTGGC" "CAGTCTCCCGCGGGGC" "GGGCCCGCCGAAAGGC"
#> [5] "GCGCGGGGCCGGTCAG" "CGGGCGCCCTGGACAC" "ACTCTGCGGGGGAGTC"

```

Created on 2023-10-12 with reprex v2.0.2

And we can doublecheck the results:

```
M <- 7000
I <- 500

[...]

# check
gc_prop <- function(x){
  toupper(x) |>
    strsplit("") |>
    sapply(\(.char) mean(.char == "C" | .char == "G"))
}

sapply(sequences, nchar) |> table()
#>
#>  500
#> 7000

min(gc_prop(sequences))
#> [1] 0.702

```

Created on 2023-10-12 with reprex v2.0.2
> 3 votes

---
Tags: r, for-loop
---


title: MegaX Pairwise Matrix Error

text:
# Question
Title: MegaX Pairwise Matrix Error
I have a pairwise matrix that I've manually formatted into a ".mega" file following the guide in the MegaX documentation but when I go to import it into the application it gives this error:

```
Unable to find the command terminator '
' (in line 48)

```

I don't understand what the terminator is supposed to be I can't find anything in the documentation or online. This is what my file looks like:

```
#mega
TITLE: Lower-left triangular matrix

#sample1
#sample2
#sample3
#sample4
#sample5
#sample6
#sample7
#sample8
#sample9
#sample10
#sample11
#sample12
#sample13
#sample14
#sample15
#sample16
#sample17
#sample18
#sample19
#sample20
#sample21
#sample22

6779.024
4010.596    5085.012
4444.925    3965.719    3597.515
3824.67     4435.878    2140.942    2855.545
2698.035    5877.724    4274.36     3721.735    3562.716
5746.71     4991.678    3103.77     4587.394    3318.78     6141.058
2898.099    7345.178    5461.126    5000.322    4846.293    2354.302    7366.433
4667.862    2949.764    3986.194    2832.598    3118.544    3551.237    5007.972    5056.938
3956.69     5110.169    2489.003    3530.347    2310.696    4285.076    3286.188    5301.894    3994.313
4086.841    5636.588    2851.673    3681.433    3078.399    4197.933    4416.542    5139.341    4347.211    3250.879
3595.868    6237.556    3924.761    3907.091    3624.064    3081.475    5770.594    3985.626    4223.117    4013.377    3136.24
4574.567    5763.568    2420.15     4433.579    3203.997    5015.929    3445.061    6063.541    4820.971    3131.669    2940.224    4132.314
4621.723    3265.735    3518.476    2528.629    2685.442    3581.311    4361.296    5144.475    2322.969    3536.469    3942.852    4218.651    4431.758
4914.659    3750.498    2719.845    3117.886    2456.598    4411.935    3294.878    5908.595    3261.36     3157.398    3657.63     4540.83     3583.543    2199.326
5662.093    3090.623    3190.918    3289.562    2772.944    5317.999    2559.264    6786.94     3426.171    3247.437    4233.766    5294.141    3960.851    2656.578    2155.925
5559.707    3073.796    4997.709    3838.601    4294.455    4486.788    5728.029    5962.611    2303.294    5051.106    5270.304    5020.468    5639.198    3504.092    4111.158    4282.95
3906.494    4028.613    3618.28     2618.348    2835.36     2743.097    4941.507    4259.628    2191.873    3599.004    3847.302    3506.872    4432.612    2061.23     2947.972    3583.198    3342.373
6582.561    1663.359    4830.862    3730.956    4161.803    5594.268    4860.298    7160.793    2675.27     4859.183    5448.895    5895.575    5537.159    3011.78     3374.505    2881.87     2979.42     3464.931
4706.715    4004.427    4618.035    3537.473    3842.462    3773.964    5703.134    5002.92     2263.311    4592.613    4872.922    4301.443    5290.578    3530.202    4175.393    4494.308    2427.841    2677.65     3599.691
5703.386    2214.959    4515.044    3273.191    3766.439    4681.361    5022.943    6197.46     1837.686    4544.365    5026.108    5096.364    5257.581    2833.977    3419.825    3267.254    2034.999    2778.249    1661.474    2332.169
6689.288    11865.223   8990.382    9113.452    8897.297    7007.562    10983.956   6142.025    9502.64     9014.072    8056.924    6431.134    8821.925    9843.732    10166.579   11038.532   10009.974   8748.955    11573.491   8952.942    10488.639

```

**EDIT:** I copied the contents of the `.meg` file over to an Excel sheet to make it easier to see but I changed the `#Mega` to `#mega` and I made sure to only include the correct elements (there are 22 samples in total, only 21 values present in the first line of the matrix). I also made sure to keep a line between the last sample name and the first value of the matrix. However, it's *still* giving me the same error.

Is there something I'm supposed to add at the end of the file to make it read the file correctly? I'm trying to import this distance matrix so I can generate a NJ tree. I'm not sure what I'm doing differently than what is mentioned here.

# Answer
Criteria, matrix. There must be a blank line between the last sample name and the start of the matrix. Also better to use #mega.

The tree from the matrix in the question is.

It's simply the carriage returns **must** by in Windows format `\r` (even if its Linux or Mac) and **not** Mac or Linux `\n`. In Mac the converter is BBEdit.
> 1 votes

---
Tags: megax, distance-matrix, pairwise, neighbor-joining
---