# Question
Title: How to calculate Allele frequency from vcf file
Organism under investigation is *Plasmodium falciparum*. How to calculate the allele frequency for each row?

I tried with this code:

```
    import pandas as pd

# Read CSV file
# df = pd.read_csv("Pf_M76611.pf7.200_samples - Copy.csv")

df= pd.read_csv("/content/Pf_M76611.pf7.200_samples - Copy.csv")

# Define function to calculate allele frequencies
def calc_allele_freq(row):
    # Get total number of samples
    num_samples = len(row) - 9

    # Get reference allele
    ref_allele = row['REF']

    # Get alternative alleles
    alt_alleles = row['ALT'].split(',')

    # Initialize list to store allele frequencies
    allele_freqs = []

    # Loop over alternative alleles
    for alt_allele in alt_alleles:
        # Initialize count for this allele
        count = 0

        # Loop over samples
        for i in range(9, len(row)):
            # Get genotype for this sample
            gt = row[i].split(':')[0]

            # Count this allele if it is present in the genotype
            if alt_allele in gt:
                count += 1

        # Calculate frequency for this allele
        freq = count / num_samples

        # Add frequency to list
        allele_freqs.append(freq)

    # Add frequency for reference allele
    allele_freqs.insert(0, 1 - sum(allele_freqs))

    # Return allele list and frequency list as a tuple
    allele_list = [ref_allele] + alt_alleles
    freq_list = [round(freq, 1) for freq in allele_freqs]
    return (allele_list, freq_list)

# Apply function to each row of the dataframe
allele_freq_df = df.apply(calc_allele_freq, axis=1, result_type='expand')

# Add columns to dataframe
df['Allele List'] = allele_freq_df[0]
df['Proportion List'] = allele_freq_df[1]

# Select columns for output and drop duplicates
# output_df = df[['#CHROM', 'POS', 'Allele List', 'Proportion List']].drop_duplicates()

# # Convert output dataframe to string and print
# print(output_df.to_string(index=False))

```

# Answer
> 0 votes
Coding: first and foremost all operations should be performed directly in `pandas`. It is much more efficient and alot easier. You don't need all the loops, `pandas` does that. `apply` is okay. I would use a very different strategy (below).

---

The actually contents of the `csv` are not clear.

This is what I think might have happened: there's multiple VCFs against the same reference I assume these were loaded and a database `join` operation performed using the reference genome column 'REF' as the index for the `join` across all input VCFs. I've placed this here because I'm not 100% sure what the input file represents - if it's not what I think it is - thats how I'd get to the allelic frequency.

Why do I think that?

1. Without multiple strains the allelic frequency cannot be calculated.
2. The file states '200 samples'
3. The presence of multiple alleles in the 'ALT' column.

The problems are:

* heterozygosity. This assumes falciparum is in its haploid state (it has a diploid stage). If its the diploid stage the allelic frequency cannot be calculated using this input file because the linkage was discarded (which is really important)
* The underlying problem is the number of actual strains isn't explicitly stated .
* the frequency of the 'ALT' doesn't appear to be declared
* even if the a frequency was declared (but not presented) that must be linked to the 'ALT' column, each allele needs a exact frequency, otherwise population genetics will not work on this data.

**First conclusion** The allelic frequency is performed across each row using the reference genome as the reference via `value_counts()` divided by the number of strains and presented in its own column. Somehow the relevant columns are missing - will need to see them.

**However** Again on the information presented I don't think it can be done because if this represents 200 samples the frequency of a given allele is not declared in this data, or simply the number of 'ALT's across all VCFs. Without accessing that number the allelic frequency is not doable because it's entirely dependent on it.

**Possible alternative conclusion** Its possible the question needs to be redefined not as *allelic frequency* \- this has an exact definition within population genetics, but *genetic diversity*, i.e. remove the term 'frequency' altogether. What I think is being calculated is genetic diversity of all samples in the file against a reference genome. If that is correct it would be good to change the question accordingly. These are different mathematical entities.

Whatever operation is required is easy to calculate, but its obtaining that information from the input file and the precise nature of the investigation.

---
Tags: python, sequence-analysis, genomics, allele-frequency
---


# Question
Title: Good / recommended way to archive fastq and bam files?
We have a lot of Illumina sequenced exome data. Currently we are using spring for its great lossless compression, but we are looking if there is anything better (and most preferably opensource) which can let us compress our fastq files.

We also want to compress or reduce file size of BAM, but using CRAM in lossless mode didn't yield good results. The space savings were decent, roughly ~10% (IMO was expecting more), but the biggest hit came when using it with IGV. We have set up a http-server and stream bams to IGV (Helped us immensely as people wouldn't have to download and visualize bam files). With Cram, the loading times were significantly high, requiring users to commit more memory for IGV. Sometimes we would have to wait around a minute for reads to show up. With Bam it is near instantaneous.

Can't use lossy compression as these BAM files might be subject to future downstream analysis (Ex. CNV calling, variant calling or something else).

What are some good ways to compress fastq and bam files so we can maximize our storage?

# Answer
EDIT: I am rewriting the answer in response to updates to the original question.

**TL;DR**: use CRAM

### Background 1: quality binning and FASTQ compression

In the old days, base callers outputted base quality at full resolution – you could see quality from Q2 to Q40 in full range. As a result, quality strings were like semi-random strings and very difficult to compress. Later people gradually realized that keeping base quality in low resolution wouldn't affect downstream analysis. The Illumina basecaller started to output quality in 8 distinct values and later changed that to 4 bins. This change greatly simplified quality string and made them compressed better. For example, in old days, a 30X BAM would take ~100 GB. With quality binning, it would only take ~60 GB.

### Background 2: GATK Base quality recalibration

In early 2010s, Illumina base quality was not calibrated well. GATK people introduced BQSR to correct that and observed noticeable improvement in SNP accuracy. Nonetheless, with improved Illumina base caller, their base quality became more accurate. Meanwhile, the world moved to 30X deep sequencing. The depth overwhelms slight inaccuracy in quality. I would say around 2015, BQSR was already unnecessary for data produced at the time.

Does it hurt to apply BQSR? Yes. First, BQSR introduces subtle biases towards the reference and towards known SNPs. Second, BQSR distorts the data. At least for some datasets, I observed that SNP accuracy dropped with variant quality after BQSR; I didn't observe this with raw quality. Third, BQSR is slow. Fourth, for new sequencers producing data at higher quality, BQSR is likely to decrease data quality. Last, related to the question, BQSR added another semi-random quality string and made compression even harder. Nowadays, running BQSR is a waste of resource for worse results. The official GATK best practice no longer uses BQSR according to their WDL file.

### CRAM

These days a 30X human CRAM only takes ~15 GB (see this file). This is a huge contrast to ~100 GB BAM in early 2010s. OP only saw ~10% saving probably due to a) BQSR and/or b) old data with full quality resolution.

On encoding/decoding speed, CRAM was much slower than BAM. Not any more. The latest htslib implementation of CRAM is faster than BAM on encoding and only slightly slower on decoding. The poor performance of IGV on CRAM could be that the Java CRAM decoder is not as optimized.

It is true that CRAM is not as widely supported as BAM. However, all the other alternatives are much worse. Petagene said they had IGV-PG (PDF) for their format. That is not official IGV and I couldn't find more recent update beyond the 2019 press release. I don't see other viable options.

Note that the common practice is to keep all raw reads, mapped or not, in BAM or CRAM such that you can get raw reads back later. BAM/CRAM additionally keeps metadata like read group, sample name, run information etc and is actually more popular than FASTQ in large sequencing centers. Also note that you don't need to sort CRAM by coordinate. Unsorted CRAM is only a little larger than sorted CRAM.

### CRAM and its competitors

The core CRAM developer, James Bonfield, is one of the most knowledgeable researchers on compression (and one of the best C programmers) in this field. He has done a lot of compression evaluation over the years. The conclusion is that on a fair benchmark, CRAM is comparable to the best tools so far in terms of compression ratio.

Petagene could compress better in the plot @terdon showed mostly because it has a special treatment of the `OQ` tag generated by GATK BQSR. It is a typical trick marketing people use to make their methods look better. With BQSR phased out, this plot is no longer relevant.

### On commercial software

In general, I welcome commercial tools and think they are invaluable to users. I also have huge respect to DRAGEN developers. However, on FASTQ storage, I would strongly recommend against closed-source compressors. If those tools go under, you may lose your data. Not worth it.
> 12 votes

# Answer
I'd recommend reconsidering CRAM. It's an open standard that is most likely to still be accessible many years in the future. Regarding your point about lossy compression...

> Can't use lossy \[CRAM\] compression as these BAM files might be subject to future downstream analysis (Ex. CNV calling, variant calling or something else).

By default, CRAM compression is not lossy.

CRAM does, however, have an optional archive settings mode (`samtools view --output-fmt cram --output-fmt-option archive ...`) which *is* a lossy compression, doing things like removing read names, removing additional accessory fields, and additional compression of quality scores. In all cases, the base sequence of the reads is preserved.

https://www.htslib.org/doc/samtools.html
> 5 votes

# Answer
You might want to consider PetaGene. This is a paid tool, and not open source, but it is something we use at work and it is really good! It can do lossless compression of gzip-compressed fastq files, bam files, and even cram files, and achieves really good size reduction (image from the PetaGene webpage):

In my own testing, using files of different sizes, I get an average size reduction of 62.1% (2.6 x smaller) for fastq.gz files (compared to the gzipped fastq file, not compared to the original, uncompressed fastq file) and an average reduction of 82.1% (5.6 x smaller) for bam files. I didn't test cram files.

If you have the budget for a paid service, these guys do good work! The greatest advantage, apart from the very significant size reduction, is that they also provide a system that allows you to use the compressed files natively. So you can feed the petagene-compressed bam file directly to bwa for example. They create special symlinks (`foo.fastq.gz -> foo.fastq.petacompressed`) and you can then use these symlinks and petagene handles the decompression on the fly.

---

The Wikipedia page on CRAM has a nice table with a few alternatives to CRAM which I reproduce below. I haven't tried any except PetaGene though so cannot comment on their performance:

| Name | Key advantage over CRAM | Repo |
| --- | --- | --- |
| Deez | None. An early approach to SAM compression, introduced 2014. | https://github.com/sfu-compbio/deez |
| Genozip | Significantly better compression, but slower random access | https://github.com/divonlan/genozip |
| GenomSys | Utilizes ISO-standard MPEG-G, but compression inferior to CRAM 3.1 | Source code not publicly available |
| PetaSuite | Virtual filesystem interface | Source code not publicly available |
> 3 votes

# Answer
An alternative to the good option of PetaGene that @terdon mentioned and also not open source is the DRAGEN ORA compression. This is loseless reference based compression and can be ~5 x smaller over fastq.gz, but it's only for fastq files and only those generated by Illumina machines (and will be easiest within the DRAGEN ecosystem).

Another point to consider is where you store the data. Even with 'standard' compression you might get more benefits moving archival data to a system like AWS S3 Glacier, depending on your current infastructure setup of course
> 3 votes

# Answer
The only free and open source tool I know that can help is `zstd`. Their github repository's README describes it as:

> Zstandard, or zstd as short version, is a fast lossless compression algorithm, targeting real-time compression scenarios at zlib-level and better compression ratios. It's backed by a very fast entropy stage, provided by Huff0 and FSE library.

We read this blog post on using it to compress fastq files and were intrigued so we ran some tests. We couldn't reproduce the same level of compression that the article claims but that might be due to different input files or the way an index was/wasn't created.

Even in our tests though, we did see that it provides better compression than gzip. Here are the specific values for one of the tests, run using an R1 file of a human WGS run:

| Uncompressed fastq | bgzip-compressed | zstd-compressed |
| --- | --- | --- |
| 179G | 58G | 47G |

This was done using `-z 19`, to set the compression level to 19 (we found this to be the best tradeoff between size and time to compress). We also saw that using a human genome fasta file to create a compression index, and then using that when compressing human fastq files, actually resulted in worse performance (larger compressed file) than not using the index. But it is entirely possible that we just did it wrong, we didn't spend a lot of time on it.

So you might want to investigate `zstd` as well, maybe you can get better performance if you play around with the indexing. On the plus side, it is free and open source. On the other hand, it won't do much with bam or cram files.
> 3 votes

---
Tags: sam, fastq, cram, storage
---


# Question
Title: compare fasta sequences in pairs and collect metrics
I have 96 fasta files (A1, A2, A3...) from one plasmid assembly pipeline, and I have another 96 fasta files (B1, B2, B3 ...) from another plasmid assembly pipeline. I would like to compare pair wise (A1 to B1, A2 to B2, A3 to B3 etc) and collect some metrics to see how this two assembly method differs.

I am thinking of doing blastn -query A1 -subject B1, but this created a lot visual alignment, rather than detailed align metrics.

Thank you for your suggestions.

# Answer
> 3 votes
It's easy; just use output `format 6`, this gives you a nice table that's easy to parse:

```
blastn -db plasmid2.db -query plasmid1.fasta -evalue 1e-6 -outfmt "6 pident evalue bitscore"

```

I don't know what parameters you prefer.

```
... "6 pident evalue bitscore" ...

```

You just edit the `pident`, `evalue` and `bitscore` above. You might want to place `staxid` or something to tell you what it is. A selection of parameters are as follows:

> staxid ssciname qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qcovs sstrand

You might want to make the command a bit fancier

```
blastn -db plasmid2.db -query plasmid1.fasta -evalue 1e-7 -outfmt "6 pident evalue bitscore" | sort | uniq -c | sort -n >> myplasmid.txt

```

---
Tags: phylogenetics, fasta, phylogeny, assembly, blast
---


# Question
Title: Parsimony-distance neighborhoods in tree space
I have the following question about the paper A parsimony-based metric for phylogenetic trees by V. Moulton and T. Wu. In this paper, the authors define a "parsimony-distance" between phylogenetic trees on the same set of taxa. They study how this distance is related to other common distances used in phylogenetics, such as counting tree bisection and reconnection (TBR) moves.

For each notion of distance, we can define the "neighborhood" of a tree as the set of trees which have distance one from the first specified tree. It is claimed in the paper that the parsimony-distance neighborhood contains the TBR-distance neighborhood, for any tree.

Lemma 6.2 states that only caterpillar trees have unit distance to some other tree. This seems to suggest that for any non-caterpillar tree, the parsimony-neighborhood should be empty. However, the TBR-neighborhood of any tree should be non-empty, since it is possible to apply a TBR move and obtain a different tree. Are there some assumptions I am missing about when the parsimony-neighborhood is empty, or when the TBR-neighborhood is contained in the parsimony neighborhood?

Also, when applying the formula in Corollary 6.6 to the non-caterpillar tree on 6 leaves, I am unsure how to apply the formula to give the neighborhood size $|N\_p| = 0$, which should be implied by Lemma 6.2.

# Answer
> 0 votes
The divergence in language makes it difficult to conceptualise and without that it's difficult to approach algebraic abstraction, albeit I've worked on tree space and its landscapes. I suspect I get point 3, but not much more.

1. "Caterpillar tree" is a term in graph theory. If that is the definition used here then every parsimony tree would be a caterpillar tree, there would never be a "non-caterpillar tree". This is because each branch length defines one "step": "parsimony's version of a mutation", so each "step" within the branch could be considered to be a "branch-less node", or else bifurcating node.

Ladderized tree I don't understand that. Maybe Np = 0 means a complete polytomy, i.e. a star tree, but could be a branch length = 0, i.e. genetically identical taxa.

2. The concept for an empty neighborhood ... possibly if all topologies are equally parsimonious this could degenerate into an empty tree neighborhood, again a star tree. What parsimony does is always define a bifurcating tree regardless, over replication all alternative solutions may result collectively in an "empty neighborhood". This might be fine from an algebraic perspective, its weird from algorithm perspective though.
3. When a TBR-neighborhood is equal to a parsimony neighborhood that actually makes sense. I define that as the TBR-neighborhood being "trapped" within the parsimony neighborhood and thats convergence. Here the "step length" has reached its minimum plateau and that tree is the most parsimonious tree ... that makes complete sense. This also makes sense within the context of the paper (I haven't read it - but I can see where they are heading).

What the authors want to do is pick holes in TBR so they want to define the point at which TBR has reached its most parsimonious tree ... I assume they want to show TBR can make errors and their solution is better. Thus they need to define the "end point". At this "end point" they want to show their measure outperforms TBR to reach a lower step size.

Thus they will want to imply that TBR should be replaced with their derivation. Its a nice methodological way to basic (or "fundamental") phylogentics.

---
Tags: phylogenetics, phylogeny
---


title: Pairwise alignment using CIGAR String

text:
# Question
Title: Pairwise alignment using CIGAR String
I want to perform pairwise alignment and get the resulting CIGAR string (Concise Idiosyncratic Gapped Alignment Report). A fuller description being here

The channel can emit insertions, deletions and substitution errors.

For example, if the input was *AAGCT*

Then the output could be *AAAAGGC*

Where there were the following errors in the following positions.

1, Insertion of A at first position 2, Insertion of A at second position 3, Insertion of G at the third position 4, Deletion of T at the last position

I want to know how I could obtain the CIGAR string given two sequences. In addition, could I use viterbi algorithm for this?

# Answer
> 1 votes
Following the basic rules and preferring the `AAG` contiguous we can write:

```
  01234
  AAGCT     (ref)
AAAAGGC     (query)

```

So the CIGAR string is `2I3M1I1D` that is:

* `2I` insert 2`A`s at the start:
* `3M` match the 3 bases `AAG`
* `1I` insert the `G`
* `1D` delete the final `T`

---
Tags: sequence-alignment, phylogenetics, errors
---


title: Singularity & Picard

text:
# Question
Title: Singularity & Picard
I am trying to run Picard in Singularity and try as I might I cannot seem to find the path to use. For example with GATK 4.4.0.0 I use the following path:

```
singularity exec $singularity_image /gatk/gatk --java-options -Xms8000m \

```

where `singularity_image` is a variable that encodes the path to the gatk.sif. However, I cannot seem to figure out the path I need to use for picard.

When I look inside the Singularity file with the following code:

```
singularity exec picard_latest.sif /bin/bash -c "ls /"

```

I get the following:

```
bin  boot  dev  environment  etc  home  lib  lib32  lib64  libx32  media  mnt  opt  proc  root  run  sbin  singularity  srv  sys  tmp  usr  var

```

**Attempt 1: Using Find**

```
singularity exec picard_latest.sif find / -name "picard.jar" 2>/dev/null

```

However, this gave no results.

**Attempt 2: Going Down The Opt Path**

```
singularity exec picard_latest.sif /bin/bash -c "ls /opt/"

```

This led to a file called `java` and that had only one file in it called `openjdk` so that seemed to be where java was stored.

**Attempt 3: Looking In Boot**

```
singularity exec picard_latest.sif /bin/bash -c "ls /boot/"

```

This yielded no results as it seems to be an empty file.

I have made other attempts but I cannot seem to establish the path I should use. I cannot change the java directly as it is preset on my cluster.

Does anyone know what the path would be for the latest version of Picard or could point me to how to look for this as I cannot seem to find this?

# Answer
> 2 votes
Assuming you've pulled the latest broadinstitute/picard container, your first attempt should have produced the expected result:

```
$ singularity exec picard_latest.sif find / -name "picard.jar" 2>/dev/null
/usr/picard/picard.jar

```

And then to run Picard and display the usage, you would use:

```
$ singularity exec picard_latest.sif java -jar /usr/picard/picard.jar -h
USAGE: PicardCommandLine <program name> [-h]

...

```

Note that there's now a newer release of Picard, which uses a much newer version of Java. There's no tag for this on Docker Hub (yet). So if you wanted the 'latest', you could try using the quay.io biocontainer:

```
$ singularity pull docker://quay.io/biocontainers/picard:3.0.0--hdfd78af_1
$ singularity exec picard_3.0.0--hdfd78af_1.sif picard -h
USAGE: PicardCommandLine <program name> [-h]

...

```

Note that the Bioconda build system also builds Docker containers and these are uploaded to quay.io as part of the Biocontainers project.

---
Tags: picard, singularity
---


# Question
Title: No common application protocol between the client and server? An error when using Entrez Direct
After a long time not running EDirect, I got an error when I tried using it. It said," curl: schannel: next InitializationSecurityContext failed: Unknown error (0×80090367) - There is no common application protocol between the client and server... I tried using WiFi and hotspot, but it did not work. Additional info: 1- I ran EDirect on Cygwin64 Terminal on my personal computer with Windows 11 operating system. 2- I first ran 'export PATH=${PATH}:${HOME}/edirect'. After that, I tried 'esearch -db pubmed -query "stroke AND english \[FILT\] AND has abstract \[FILT\] NOT review \[PTYP\]" | efilter -mindate 2000 | efetch -format xml | xtract -pattern PubmedArticle -element PMID AbstractText \> dm\_cs.txt'

I used to ran it successfully on the same PC over WiFi or hotspot three month ago.

# Answer
This happens.

1. One possibility is it's simply NCBI blocking the request to avoid being overrun by e.g. bots. This occurs if you've been running the script repeatedly within a single day. It will be fine tomorrow in this scenario.
2. The second possibility is you've not supplied an email address.
3. The final possibility is if you're running a VPN. NCBI has taken a strong dislike to VPN's *per se*, I've never figured out its algorithms - it seems to permit a fixed number of requests then pointedly refuses. Just turn off the VPN.

If its point 2 then understanding how you've coded this would be useful to allow an email address to be sent. Biopython has a function for this. I just use Biopython BTW.
> 1 votes

---
Tags: curl, edirect
---


# Question
Title: How much of the space of motifs have been discovered?
I notice when using tools like HOMER there are two different types of outputs: de novo motifs vs known motifs and it seems no matter what the data is (type of ChIP etc.) there is always a "de novo" motif set in the results.

I was hoping to reach out to experts on insight about how much of the space of motif based structure, within a genome, has been identified. Is there any consensus on what portion of the motif like patterns, for example in human or mouse genome, have been identified? When I find de novo motifs using HOMER in my ChIP-seq or ATAC peaks, what are the chances that in my common-place ChIP experiment (PU.1 for example) I find de novo motifs, not recorded in the databases (i.e. JASPER).

# Answer
> 2 votes
De novo means that the software infers presence of motifs based on the provided sets of DNA sequences, most commonly peaks from ChIP-seq and similar experiments. So essentially it stacks the sequences on top of each other and checks for presence of conserved nucleotides that form sequences of certain lengths, for example \< 10bp.

Known means that each peak is checked whether it contains a known motif (so is similar to a given known motif) and then in the end a per-motif statistic is calculated to decide whether enriched or not.

If the investigated protein is well-studied then de novo and known is likely similar or identical. De novo does not mean that found motif is novel in terms or unknown or undescribed. It just means that it is derived from the data directly rather than matching known motifs onto the peaks. It's simply two different ways to approach the question. Hope that makes sense.

---
Tags: motifs, homer
---


title: Problem with Seurat reference mapping

text:
# Question
Title: Problem with Seurat reference mapping
I have `10X scRNA-seq multiome's 3' poly-A capture (scRNA-seq+ATAC-seq)` from PBMC

Using Seuratreference mapping, I mapped my scRNA-seq part on reference PBMC (for cell type annotation)

I have also tried Seurat zimuth shiny app

But likely just a few of my cells has been mapped on the reference

Please look at this

Do you think this because my data is multiome or something else is going wrong?

Please any body has some experience help me

Thank you

# Answer
> 1 votes
Yes! Multiome seems to have cells with much higher mitochondrial gene content. I don't understand the biology behind this. But in using a cut off of 5% (i.e. remove cells with \>5% MT genes) I lose over 80% of cells in multiome data. Also, here is the Signac tutorial where they use a 20% cut off for multiome data - https://satijalab.org/seurat/articles/weighted\_nearest\_neighbor\_analysis.html ("WNN analysis of 10x Multiome, RNA + ATAC")

# Answer
> 0 votes
As I mentioned in my answers to similar questions, there are a lot of cells with very high mitochondrial gene content. It would be a good idea to work out why that is before doing any more downstream processing.

These cells will, of course, be filtered out by the thresholds you have set.

---
Tags: scrnaseq, seurat, sequence-annotation, 10x-genomics, cellranger
---


# Question
Title: calculate mismatch frequency/rate from a BAM file
I am working with PAR-CLIP data where experimental procedure induces T to C mismatches. The other types of mismatches we can consider coming from artifact such as PCR errors. It could also be an existing SNP. For simplicity, let's consider T to C mismatches are experimentally induced and rest of the mismatches are background mismatch rate.

From Illumina NGS seq after performing alignment, for each of he `read groups` (reads collapsed in a 30-70 base region with minimum 10 read coverage). Up to this point its already done.

I want to compare T to C mismatches vs background mismatches. To do so I need to calculate the rate of signal mismatches (T to C) and background mismatches.

Using `samtools mpileup` I can get total number of all kind of mismatches for a `read group`, however, takes some time.

1. How do I calculate background mismatch rate from this?
2. How do I calculate total mismatches in a bam?

I was thinking about something like
for a read group: `mismatch_rate_group=total_mismatches_group/(total_read_group*total_base_group)` for a BAM file: `mismatch_rate_global=total_mismatches_bam/(total_reads_aligned_bam*genome_size_base)`

3. What important factor I might missing in my approach?

I am not using any variant caller (e.g. bcftools). Initially I tried with samtools mpileup and then used a shell script to figure out the mismatches. But I am not sure is this is the right thing to do. Or even there are any easier methods. If one read has a mutation, I consider that a mismatch and I presume it to be artifact and/or PCR errors.

# Answer
> I am not using any variant caller

Use a variant caller. `bcftools`, at least, will output variant frequency per variant:

```
bcftools mpileup -ugf ref.fa sample.bam | bcftools call -mv > output.vcf

```

As an alternative approach, I wrote some code to parse the mpileup output and extract counts for different alleles; this output may be easier for you to parse and calculate your desired frequencies:

```
# readstomper.pl -- generate base-packed proportion statistics for
#   output from 'samtools mpileup'
#
# example use:
# samtools mpileup -B -Q 0 -f circ-Nb-ec3-mtDNA.fasta LAST_OL_132394_vs_mtDNA.bam | \
#   /bioinf/scripts/readstomper.pl > LAST_OLmpileupProp_132394_vs_mtDNA.txt.gz
#
# example output:
# Assembly,Position,Coverage,ref,cR,pR,A,C,G,T,d,i,InsMode
# Nb_mtDNA,11456,86,T,79.00,0.92,0.0000,0.0465,0.0000,0.0000,0.03,0.01,A;1.00
# Nb_mtDNA,11457,86,A,81.00,0.94,0.0000,0.0000,0.0233,0.0000,0.03,0.02,TA;0.50
# Nb_mtDNA,11458,86,G,79.00,0.92,0.0349,0.0233,0.0000,0.0116,0.01,0.00
# Nb_mtDNA,11459,86,A,78.00,0.92,0.0000,0.0000,0.0235,0.0471,0.01,0.00
# Nb_mtDNA,11460,86,T,5.00,0.06,0.0233,0.8140,0.0233,0.0000,0.08,0.00
# Nb_mtDNA,11461,86,G,57.00,0.66,0.1512,0.1163,0.0000,0.0465,0.02,0.01,AATAACAACACGTAACCGA;1.00
# Nb_mtDNA,11462,86,G,78.00,0.91,0.0349,0.0349,0.0000,0.0116,0.01,0.00
# Nb_mtDNA,11463,86,G,66.00,0.77,0.0465,0.0698,0.0000,0.0116,0.10,0.01,GCAATA;1.00
# Nb_mtDNA,11464,86,T,74.00,0.86,0.0000,0.0814,0.0233,0.0000,0.03,0.00

```
> 0 votes

---
Tags: ngs, sam, point-mutation
---


title: cmapPy exception

text:
# Question
Title: cmapPy exception
I'm trying to read gctx files with cmapPy. This is my code:

```
from cmapPy.pandasGEXpress.parse_gct import parse
from cmapPy.pandasGEXpress.write_gct import write
parse("C:/Users/User/Desktop/gctx/CAS9_BASELINE_LEVEL3_BING__n10721x10174.gctx")

```

This is the error:

```
(base) C:\Users\User\Desktop\hello>D:/anaconda/python.exe c:/Users/User/Desktop/gctx/readgcxtx.py
Traceback (most recent call last):
  File "c:/Users/User/Desktop/gctx/readgcxtx.py", line 3, in <module>
    parse("C:/Users/User/Desktop/gctx/CAS9_BASELINE_LEVEL3_BING__n10721x10174.gctx")
  File "D:\anaconda\lib\site-packages\cmapPy\pandasGEXpress\parse_gct.py", line 136, in parse
    num_row_metadata, num_col_metadata) = read_version_and_dims(file_path)
  File "D:\anaconda\lib\site-packages\cmapPy\pandasGEXpress\parse_gct.py", line 168, in  read_version_and_dims
    version = f.readline().strip().lstrip("#")
  File "D:\anaconda\lib\encodings\cp1252.py", line 23, in decode
    return codecs.charmap_decode(input,self.errors,decoding_table)[0]
UnicodeDecodeError: 'charmap' codec can't decode byte 0x90 in position 2440: character maps to <undefined>

```

I tried changing the path with double slashes or r"..." but it doesn't work. Thanks.

# Answer
> 1 votes
The error I see is,

> UnicodeDecodeError: 'charmap' codec can't decode byte 0x90 in position 2440: character maps to

This suggests there is an error in the input format that the program doesn't like.

The thing that bothers me is why you have a mix of / and \ for directory paths. Windows uses \ and Linux and OS X use / . I don't Windows so I really can't comment.

1. It looks like you've used a character at position 2440 that the code doesn't like, this is often the case with ambiguous characters such as degenerate nucleotides or ambiguity codes for example "?" or "N". It seems to suggest its neither or these but something like ...
2. I'd be cautious of is the carriage returns because Windows is of course very different to Linux in this regard. \r (I think) versus \n

`undefined` mapping suggests is quite a specific error and therefore resolvable.

# Answer
> 0 votes
from cmapPy.pandasGEXpress.parse import parse

from cmapPy.pandasGEXpress.write import write

parse("C:/Users/User/Desktop/gctx/CAS9\_BASELINE\_LEVEL3\_BING\__n10721x10174.gctx")

---
Tags: python, anaconda
---


title: Strategy for merging many VCFs

text:
# Question
Title: Strategy for merging many VCFs
Can anyone suggest a strategy for speeding up VCF merging?

I have ~44,000 single-sample VCFs that I am trying to merge into a multi-sample VCF with `bcftools merge`, but the job keeps timing out on the cluster (just timed out after 4 days). I could obviously try running it for a couple of weeks, but I was hoping there is a better strategy someone can suggest?

The exact command I am running is

```
bcftools merge --no-version -m all -i 'DP:avg,RO:avg,AO:avg' -F + -o merged.vcf --file-list merge.fofn

```

Happy to try any other tool that will give the same merge results.

# Answer
We created this over the pandemic when we needed to merge loads of VCFs which all have the same header and first few columns (ie same variants)

https://github.com/iqbal-lab-org/ivcfmerge

That page gives detailed instructions but i'll paste a couple of examples of how you can use it

== In Python =4.1.1 If the number of input files is small (can be opened all at once)

```
from contextlib import ExitStack
from ivcfmerge import ivcfmerge

filenames = [...]    # List/iterator of relative/absolute paths to input files
output_path = '...'  # Where to write the merged VCF to

with ExitStack() as stack:
    files = map(lambda fname: stack.enter_context(open(fname)), filenames)
    with open(output_path) as outfile:
        ivcfmerge(files, outfile)

```

=4.1.2 If the number of input files is big (cannot be opened all at once) from ivcfmerge import ivcfmerge\_batch

```
filenames = [...]    # List/iterator of relative/absolute paths to input files
output_path = '...'  # Where to write the merged VCF to
batch_size = 1000    # How many files to open and merge at once

ivcfmerge_batch(filenames, output_path, batch_size)

```

Pretty sure we initially did batched bcftools merge, I can't honestly remember why we moved to this. Thankfully not having to do these merges very often
> 4 votes

# Answer
I ended up using an iterative approach (inspired by this post), which took 1.5 hours - so a massive speed up (and that was just the sequential version).

Note: `fd` is analogous to `find`, but way nicer to use.

First, merge subsets of *n* VCFs

```
# size of the subsets
n=200
# create file of filenames for each subset
fd -e vcf . path/to/vcfs/ | split -l $n --additional-suffix '.fofn' - subset_vcfs

```

`-e vcf` tells `fd` to find all files with a `.vcf` extension, you can obviously change this to `.bcf`/`.vcf.gz` etc. if necessary. `split` splits the list of all VCFs files into subsets of size `-l` (200) and creates the file of filenames with the prefix `subset_vcfs`. It will add a random suffix to this - e.g. `subset_vcfsav`. Plus the additional suffix `.fofn`. See the `split` manual for all options.

Now iterate through these subsets and merge the VCFs listed in them.

## Subset merge option 1 - sequential

This just merges the VCFs one at a time

```
for sub in *.fofn;
do
    echo "Merging subset ${sub}..."
    out_vcf="merge.${sub}.bcf"
    bcftools merge --file-list "$sub" -o "$out_vcf"
    bcftools index "$out_vcf"
done

```

## Subset merge option 2 - parallel

You can easily parallelise the previous step using `fd` on the subset files. See this post for more details. Or basic bash `for` loop parallelisation.

Create a script (`merge_subset.sh`) of what we want to run on each subset

```
#!/usr/bin/env bash
sub="$1"
echo "Merging subset ${sub}..."
out_vcf="merge.${sub}.bcf"
bcftools merge --file-list "$sub" -o "$out_vcf"
bcftools index "$out_vcf"

```

Now we run this script with each subset file of filenames using `fd`

```
# number of processes to use. remove this option to use all available
jobs=4
fd -j $jobs -e fofn -x bash merge_subset.sh '{}'

```

## Final merge

Now we merge all of these VCFs into our final VCF

```
ls merge.*.bcf > merged.fofn
bcftools merge -o merged.bcf --file-list merged.fofn

```
> 3 votes

# Answer
Perform a tree based merge of the files.

* Sort the files into groups that can easily be merged in a measurable amount of time.
* Merge each group into a separate Variant Call Format (VCF) file
    + Has the benefit that files can also be grouped by type / source / run
    + Further benefit that small merges can be checked to ensure you like the output before proceeding to the next stage.
    + Further benefit that there are also subgroup files available to be distributed or shared separately.
    + Downside, extra storage space until you have a final product.
* Merge subfiles upward in stages until you arrive at a complete 44,000 entry file
    + If you're happy with the earlier test results, process can be automated with relatively simple iterative shell scripts.
    + Possibly faster, as the bulk merge needs to load every single file into memory, and you may be hitting a capacity bottleneck.
    + Worked in Fluid Dynamics for years, and often ran into a similar issue. Could run many small cases representing a large case faster than running the large case, because large case would overwhelm memory and transfer capacity.
    + In some clusters you may also be causing nodes to silently crash, which then stalls the entire process without a visible notification.
> 1 votes

---
Tags: vcf, bcftools, merge
---


title: Manta Evidence Folder Is Empty

text:
# Question
Title: Manta Evidence Folder Is Empty
I am looking at Manta on Github and I ran a bit of code:

The following variables were defined in the code:

* `input_dir` is the file path to a folder with less than 10 .bams as per the limitations of Manta
* `output_dir` is the path to the output
* `MANTA_INSTALL_PATH` is the path to the Manta installation
* `reference` is the path to the reference fasta

```
    for file in ${input_dir}/*.bam
do
  # Create subdirectory named after bam file
  sub_dir=${output_dir}/$(basename ${file} .bam)
      mkdir -p ${sub_dir}

      # Use reference fasta path specified by user
      ${MANTA_INSTALL_PATH}/bin/configManta.py \
    --bam $file \
        --referenceFasta ${reference} \
    --runDir ${sub_dir} \
        --exome
    done

# Now run each `runWorkflow.py`

find "$output_dir" -name "runWorkflow.py" -exec {} \;

```

The run worked and each sample had a folder with the following files:

```
results  runWorkflow.py  runWorkflow.py.config.pickle  workflow.error.log.txt  workflow.exitcode.txt  workflow.warning.log.txt  workspace

```

Each `workflow.error.log.txt` had nothing in it implying the the run had been successful. Each `results` folder had 3 subfolders

```
evidence  stats  variants

```

The `stats` and `variants` folder each had outputs but the `evidence` folder had nothing. Is this normal? Should it have something in this folder? Why is this folder there in Manta if nothing is in it?

# Answer
> 1 votes
The documentation on the GitHub page you linked to has this to say about the `evidence` folder:

> **Generating evidence bams**
>
> Using the `--generateEvidenceBam` option, Manta can be configured to generate bam files of evidence reads for SVs listed in the candidate vcf file.
>
> It is recommended to use this option together with the --region option, so that the analysis is limited to relatively small genomic regions for debugging purposes.
>
> The evidence bam files are provided in `${MANTA_ANALYSIS_PATH}/results/evidence`, with a naming format `evidence_*.*.bam`. There is one such file for each input bam of the analysis, containing evidence reads of the candidate SVs identified from that input bam. Each read in an evidence bam keeps all information from the original bam, and it contains also a customized tag in the format: `ZM:Z:${MANTA_SV_ID_1}|${EVIDENCE_TYPE},${MANTA_SV_ID_2}|${EVIDENCE_TYPE}`. For example, `ZM:Z:MantaINV:5:0:1:0:0:0|PR|SRM,MantaDEL:5:1:2:0:0:0|SR`
>
> One read can have more than one of the three evidence types: PR for paired reads, SR for split reads, and SRM for split read mates. One read can be evidence for multiple SVs, which are separated by commas in the tag. Notice that the number of evidence reads for a particular SV in the evidence bam files could be more than the evidence counts (PR and SR) in the final vcf files. This is because more stringent criteria are applied for generating evidence counts in the final vcf files.

So it looks to me like you need the `--generateEvidenceBam` option to generate evidence, but that the directory `evidence` is created regardless of whether you use the option. Try running the analysis again, but using `--generateEvidenceBam` this time.

---
Tags: git
---


# Question
Title: Integrating bulk RNA-Seq data with different sequencing depths and from different sources
I am attempting to integrate different bulk RNA-Seq datasets. While this is not ideal, I'm trying to reduce the technical variability in these datasets by using data generated by similar protocols (for example, all of them have been sequenced with a Poly-A enrichment strategy and Illumina short read sequencing). I am aware the results I obtain from this analysis will be very limited, but still want to pursue the analysis. I'm looking for advice on two aspects of the analysis, but first some more info about the data:

* There are 5 bulk RNA-Seq datasets and 4 biological groups. Each datasets only comprises one biological group. Only two biological groups have two datasets, and such datasets overlap in PCA.
* In Principal Component Analysis (PCA) and the first PC (~26% of variation) is very correlated with the biological variable (cell differentiation), suggesting there is a biological signal among all the noise/technical bias.

Despite these preliminary results I'm not sure I'm accounting correctly for some aspects: The datasets have been sequenced with different sequencing coverage (in fact, different numbers and length of reads). Some are single read and some paired-end. My questions are:

1. For PCA, I'm normalizing with variance stabilizing transformation (VST). How should I account for the differences in read numbers and type between datasets?
2. For differential expression (DE) analysis I'm using DESeq2. How should I account for the differences in read numbers and type between datasets? There is no batch variable, as all datasets come from different sources (batches, laboratories)

I'm happy to provide more information if necessary. Thanks

# Answer
Per-sample coverage / count differences are expected in all differential expression analysis, and correction for that is built into DESeq2. You shouldn't need to do any additional correction with respect to proportional counts when looking at a PCA derived from the VST.

For accounting for different sources / types / etc, it's a good idea to add that information into the metadata, and test for differential expression with and without accounting for those variables to see if anything changes. This is recommended in the DESeq2 manual:

> If there is unwanted variation present in the data (e.g. batch effects) it is always recommend to correct for this, which can be accommodated in DESeq2 by including in the design any known batch variables or by using functions/packages such as svaseq in sva (Leek 2014) or the RUV functions in RUVSeq (Risso et al. 2014) to estimate variables that capture the unwanted variation. In addition, the ashr developers have a specific method for accounting for unwanted variation in combination with ashr (Gerard and Stephens 2017).

See the multi-factor design section of the manual for more information (which has a "single-read" vs "paired-end" example):

https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#multi-factor-designs

If after accounting for the different variables there is not enough information to properly separate groups, DESeq should tell you. Ideally you should aim for six replicates of each state of each independent variable (easier to do when there are interacting variables), but it's quite rare to see that in differential expression analysis.
> 2 votes

# Answer
In case of different sequencing depth one performs *rarefaction*, to make the data comparable.

The standard approach is *rarefaction by size*, when one takes the number of reads in the smallest sample, and samples without replacement the same number of reads from the other samples. Potentially, one could carry out upsizing to a greater number of reads (e.g., if the smallest sample is an outlier in terms of the number of reads, and downsizing may result in losing too much information.) I suggest reading the tutorial for MOMR package for more detailed description (the R code is also available.)

There are more sophisticated methods of rarefaction, such as *rarefaction by coverage* \- these are of significant scientific interests, but tricky to use and not very widespread. See, e.g., this article (it deals with sampling species, but generalization to reads is straightforward... though reads are not species.)

**Remark:** Of course many good packages (like those mentioned by @gringer) might already perform rarefaction by default - but it is always a good idea to check the documentation.
> 0 votes

---
Tags: rna-seq, rna, quality-control, data-preprocessing, public-dataset
---


# Question
Title: I am trying to find correlation between two ChIP-seq datasets but the Pearson and Spearman test results do not correspond with the peak overlap
I am trying to compare two ChIP-seq datasets for histone modifications, however, the Spearman and Pearson correlations between two particular datasets does not correspond to their peak overlap (peak overlap is very low, while both pearson and spearman tests are showing strong correlation). I am performing this analysis on Usegalaxy, first using multibamSummary on both bam files followed by plotCorrelation. What could be the error here?

# Answer
Correlations like this between ChIP-seq datasets is common.

If you're calculating the correlation between BAM files across the entire genome, those are largely going to be the same (hence the need for input controls in ChIP-seq for performing the peak calling in the first place). Similarly, many ChIP-seq experiments only have 30 - 60% of reads within their peaks (FRiP), which means that 40 - 70% of reads don't fall within peak regions. This is going to heavily skew any correlation calculations that don't exclusively focus on the peaks themselves.

To more clearly compare ChIP-seq data against each other, global metrics like FRiP can help you diagnose measurement/sequencing issues. Looking at the correlation across ChIP-seq signal within peaks may be a more informative calculation than looking at all reads across the genome. Similarly, principal component analysis on signal within peaks can be helpful here.

DiffBind has some helper functions to perform these types of calculations that you can look at.
> 1 votes

---
Tags: chip-seq, correlation
---


# Question
Title: Error in FindIntegrationAnchors- Seurat package
I am working on integrating a labelled single cell RNA seq cell atlas with an unlabelled one. I am wondering how do I determine the value of max.features when integrating two single cell RNA seq datasets?

Below is the code I am using:

```
epithelial = Read10X(data.dir="/projects/b1025/sdi0596/Covid19_Single_Cell_RNASeq_project/filtered_matrices_count/epithelial/filtered_feature_bc_matrix/",gene.column=1,cell.column=1)
habermann = Read10X(data.dir="/projects/b1025/sdi0596/Covid19_Single_Cell_RNASeq_project/filtered_matrices_count/Habermann_Single_Cell_Data/",gene.column=1,cell.column=1)

data.list=c(epithelial,habermann)

normalize and identify variable features for each dataset independently
ifnb.list <- lapply(X = data.list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

Perform integration
immune.achors= FindIntegrationAnchors(object.list = ifnb.list,max.features = 200,verbose = TRUE)

Error in FindIntegrationAnchors(object.list = ifnb.list, max.features = 200, :
Max dimension too large: objects 1, 2 contain fewer than 30 cells.
Please specify a maximum dimensions that is less than the number of cells in any object (4).

```

Regarding the dimension of these objects, the epithelial one is 33538 x 49109 and the habermann is 33694 x 220213. Any insights would be highly helpful.

Thank you.

# Answer
> 2 votes
There's a few problems with your code, first, when you do `Read10X()` it returns you a sparse matrix, and you need to put this into a `Seurat` object with meta data, before doing the integration.

So for example, i create some example data that is similar to your output from `Read10X()` :

library(Matrix)

```
epithelial = Matrix(matrix(rnbinom(100000,mu=200,size=1),ncol=50),sparse=TRUE)
rownames(epithelial) = paste0("gene",1:2000)
colnames(epithelial) = paste0("e",1:50)

habermann = epithelial + rnbinom(100000,mu=20,size=1)
colnames(habermann) = paste0("h",1:50)

```

Then convert to a seurat object:

```
data.list=c(epithelial = CreateSeuratObject(epithelial),
            habermann = CreateSeuratObject(habermann))

```

Then find variable features:

```
data.list <- lapply(data.list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

```

Should look like this:

```
$epithelial
An object of class Seurat
2000 features across 50 samples within 1 assay
Active assay: RNA (2000 features, 2000 variable features)

$habermann
An object of class Seurat
2000 features across 50 samples within 1 assay
Active assay: RNA (2000 features, 2000 variable features)

```

Now perform integration, below I have to reduce `k.filter` because I have very little cells in this example. In your case, you can simply use the default settings.

```
features <- SelectIntegrationFeatures(object.list = data.list)
anchors= FindIntegrationAnchors( data.list,max.features = 200,
k.filter=50,k.anchor = 3,verbose = TRUE)

```

# Answer
> 0 votes
It is possible that some of your samples have too few cells, so it is recommended to check the data and remove samples with low cell counts.

---
Tags: scrnaseq
---


# Question
Title: What does the IGV refseq genes with three different thicknesses of lines mean?
What does the IGV refseq genes with three different thicknesses of lines mean?

# Answer
> 5 votes
The thinkness indicates whether something is an intron, a coding exon, or a UTR. The very thin lines are introns, the full thickness boxes are coding exons and the thinner boxes are UTRs.

UTRs, untranslated regions, are exons, they are included in the final mRNA, but they don't code for protein. So, by convention, they are shown as thinner exons to give the user a visual indication that although this is still an exonic region, it's a special one.

I don't really know where this is documented, it's pretty standard across all genome browsers and bam visualizers I know though so it's a common idiom.

---
Tags: igv
---


# Question
Title: DESeq2 throwing error while normalizing raw microarray expression data due to presence of negative values
*This question was also asked on Biostars*

I am trying to download and analyze a miRNA expression dataset from NCBI GEO `(GSE25631)`. I specifically want non-normalized data to perform normalization and my own set of other analyses later on. Accordingly, I downloaded the `Series Matrix File` (to match the Sample IDs with their phenotype) and the non-normalized data `GSE25631_non-normalized_data.txt.gz` programmatically using the `GEOquery` package in R. To normalize the data I'm using the `DESeq2` package. When I run the `DESeqDataSetFromMatrix` function I get this error:

> Error in DESeqDataSet(se, design = design, ignoreRank) : some values in assay are negative

I went back to the non-normalized data and observed that there are a lot of negative values. I need suggestions on how to deal with this issue. Should I replace all such values with NA? I also observed that every probe has a detection p-value associated with it. Should I also filter out the expression values which have a p-value higher than 0.05?

```
> head(raw_counts)
        gbm-0001 Detection Pval gbm-0011 Detection Pval gbm-0013 Detection Pval gbm-0024 Detection Pval gbm-0058 Detection Pval
        ILMN_3167304   27.875      0.3071961 -24.9375      0.6689599    14.25       0.350027   4.3125      0.4723163     8.25      0.3966338
        ILMN_3168038  211.875       6.42E-05   1.0625      0.4925717    99.25    0.003645935  60.3125      0.1657064    72.25     0.01086376
        ILMN_3167890 13839.88       3.68E-38 23507.06       3.68E-38 22345.25       3.68E-38 20482.31       3.68E-38 21419.25       3.68E-38
        ILMN_3167526    5.875      0.4577177 105.0625     0.03279042     8.25      0.4117529 -29.6875       0.683706    11.25      0.3604082
        ILMN_3167209   82.875     0.06708142 182.0625      0.0007096    60.25     0.05167192   9.3125      0.4403947    72.25     0.01086376

```

Here's the R script which I used for all the mentioned steps.

```
library(GEOquery)
library(DESeq2)
library(tidyverse)

data <- getGEO(GEO = "GSE25631",
           destdir = "E:\\GSE25631",
           GSEMatrix = TRUE,
           AnnotGPL = FALSE,
           getGPL = FALSE,
           parseCharacteristics = TRUE)

clindata <- data$GSE25631_series_matrix.txt.gz@phenoData@data

url = "https://www.ncbi.nlm.nih.gov/geo/download/?
acc=GSE25631&format=file&file=GSE25631%5Fnon%2Dnormalized%5Fdata%2Etxt%2Egz"
download.file(url, "E:\\GSE25631\\GSE25631_raw_counts.gz")
raw_counts <-read.delim("E:\\GSE25631\\GSE25631_raw_counts.gz",
                     stringsAsFactors = FALSE, sep = "\t")
colnames(raw_counts) <- raw_counts[5,]
raw_counts <- raw_counts[-(1:5),]
raw_counts <- raw_counts %>% remove_rownames %>% column_to_rownames(var = "ID_REF")
raw_counts <- raw_counts[ , c(TRUE, FALSE)]

colnames(raw_counts) <- clindata$title
rownames(clindata) <- clindata$title
colnames(clindata)[colnames(clindata) == "disease state:ch1"] <- "Condition"
clindata$Condition[clindata$Condition == "primary glioblastoma multiform (brain tumor)"] <- "GBM"
clindata$Condition[clindata$Condition == "normal"] <- "Normal"
clindata$Condition <- as.factor(clindata$Condition)

dds <- DESeqDataSetFromMatrix(countData = raw_counts,
                          colData = clindata,
                          design = ~ Condition)

```

# Answer
> 2 votes
Answered at biostars: https://www.biostars.org/p/9561238/#9561240

> DESeq2 is meant for digital counts (RNA-seq and similar assays). Arrays are not supported, nor is the implemented method meaningful for array intensity values.
>
> Guides on how to handle these data are probably in here:

---
Tags: deseq2, normalization, microarray, geoquery, geo
---


title: Can bwameth accept unaligned Bam?

text:
# Question
Title: Can bwameth accept unaligned Bam?
Does bwameth accept `bam` (unaligned) type as input or I will have to convert unaligned bam to fastq (`samtools fastq bam`) to use it as input?

# Answer
No, it's a bwa mem wrapper for bisulfite data essentially and bwa mem accepts fastq only.
> 2 votes

---
Tags: sam, bwa, methylation
---


# Question
Title: How do i search by these gene properties in NCBI's GenBank?
I want to search GenBank entries by gene product and taxon id.

Taking the 16s rRNA as an example for a gene and Piscirickettsia salmonis (taxon:1238) for a taxon:

> https://www.ncbi.nlm.nih.gov/nuccore/219846389

How do I retrieve all GenBank entries that have the following properties **/product**="16S ribosomal RNA" and **/db\_xref**="taxon:1238"?

Update: It would be appreciated if you also showed me how to retrieve GenBank entries like this, which don't have an exact match to the "16S ribosomal RNA" product.

# Answer
This can be done using the NCBI eutils command-line tools. For taxa, the correct incantation for the search string is `txid<taxID>[Organism]` (assuming you want that taxon and all its children), but I haven't found the right one for the product field. However, it can be approximated by searching for a string in the tile line via `<search>[TITL]`. This seems to work fairly well, assuming you want the fasta files (the `-format` argument can be changed to `genbank` if you want the genbank record instead):

```
$ esearch -db nuccore -query "txid1238[Organism] AND 16s[TITL]" | efetch -format fasta | grep -A 1 '^>' | head
>MW433822.1 Piscirickettsia salmonis strain HRV16 16S ribosomal RNA gene, partial sequence
TCAGATTGAACGCTGGTGGCATGCTTAACACATGCAAGTCGAACGGTAACATGACTAGCTTGCTAGTTGA
--
>MW433821.1 Piscirickettsia salmonis strain HRV13 16S ribosomal RNA gene, partial sequence
AATTTCAGGTGTAGCGGTGAAATGCGTAGATATCTGAAGGAACACCGGTGGCGAAGGCGACTTTCTGGAT
--
>MN023100.1 Piscirickettsia salmonis clone Ps18621-6 16S ribosomal RNA gene, partial sequence
GCTTGCTAGATGCTGACGAGTGGCGGACGGGTGAGTAACGCGTGGGAATTTACCTTTTAGTGGGGGATAA
--
>MN023099.1 Piscirickettsia salmonis clone Ps18429-3 16S ribosomal RNA gene, partial sequence

```

Unfortunately I can't find anything for searching specifically in the features of an NCBI record. There'll be some way to do it; I just don't know how.

Based on the updated request, it's possible to do a full-text search through the entire metadata (e.g. for `16S ribosomal RNA`) by adding the search term without specifying a field. Depending on what you want to do with the results, and how many were produced, it could be easy enough to manually unwanted sequences (e.g. whole genome shotgun assemblies, or chimeric products):

```
$ esearch -db nuccore -query "txid1238[Organism:exp] AND (16S ribosomal RNA)" | efetch -format fasta | grep -A 1 '^>' --no-group-separator | tail -n 12
>EU170254.1 Piscirickettsia salmonis strain NVI 5692 16S ribosomal RNA gene, partial sequence
CACATGCAAGTCGAACGGTAGCAGATCTAGCTTGCTAGATGCTGACGAGTGGCGGACGGGTGAGTAACGC
>U36941.1 Piscirickettsia salmonis 16S ribosomal RNA gene, complete sequence
AGAGTTTGATCCTGGCTCAGATTGAACGCTGGTGGCATGCTTAACACATGCAAGTCGAACGGTAGCAGAT
>U55015.1 Piscirickettsia salmonis 16S ribosomal RNA gene, complete sequence
AGAGTTTGATCCTGGCTCAGATTGAACGCTGGTGGCATGCTTAACACATGCAAGTCGAACGGTAGCAGAT
>U36942.1 Piscirickettsia salmonis strain NOR-92 16S ribosomal RNA gene, complete sequence
AGAGTTTGATCCTGGCTCAGATTGAACGCTGGTGGCATGCTTAACACATGCAAGTCGAACGGTAGCAGAT
>U36940.1 Piscirickettsia salmonis strain EM-90 16S ribosomal RNA gene, complete sequence
AGAGTTTGATCATGGCTCAGATTGAACGCTGGTGGCATGCTTAACACATGCAAGTCGAACGGTAGCAGAT
>U36915.1 Piscirickettsia salmonis 16S ribosomal RNA gene, complete sequence
AGAGTTTGATCATGGCTCAGATTGAACGCTGGTGGCATGCTTAACACATGCAAGTCGAACGGTAGCAGAT

```

FWIW, here's the list of all the searchable fields in the `nuccore` database:

```
$ einfo -db nuccore -fields
ACCN    Accession
ALL     All Fields
ASSM    Assembly
AUTH    Author
BIOS    BioSample
BRD     Breed
COMP    Component Accession
CULT    Cultivar
DIV     Division
ECNO    EC/RN Number
FILT    Filter
FKEY    Feature key
GENE    Gene Name
GPRJ    BioProject
ISOL    Isolate
ISS     Issue
JOUR    Journal
KYWD    Keyword
MDAT    Modification Date
ORGN    Organism
PACC    Primary Accession
PAGE    Page Number
PDAT    Publication Date
PORG    Primary Organism
PROP    Properties
PROT    Protein Name
SLEN    Sequence Length
SQID    SeqID String
STRN    Strain
SUBS    Substance Name
TITL    Title
UID     UID
VOL     Volume
WORD    Text Word

```
> 1 votes

# Answer
To specifically answer the question the search that is used is

```
((1000:1500[Sequence Length]) AND Piscirickettsia salmonis[Organism]) AND 16S

```

That is place here and select 'nucleotide'.

Genbank uses taxon id internally, it's a database reference not a biological code, and can't be used as input - its completely unique to Genbank. Anyway the search results in 276 records, and are shown here, this is near full length 16S thats cool.

**Clarification** `[taxon id]` is not a field in Genbank, it's a database code, the tag is `[organism]`. Genbank knows what 16S is - its the most famous genetic marker in microbiology - so there is no need for a `[tag]`

At the base of the page there is a "Send to:" option (right hand corner). Click and then select "File". The 276 sequences will then download.

If `[Sequence length]` is removed this result 3820 sequences here. There are some 16S-23S sequences on the database result in a 2,198 bp length, so extending the sequence length will trap these by `1000:2200[Sequence Length]`. There are 431 bp sequences - which are not cool at all, probably a lot.

The approach used below identifies nearly 1000 is complicated it traps everything of relevance which is close to full length and cuts out the 16S from 2,198 bp length and beyond. The 431 bp partial 16S sequences have been removed. Thus the numbers below are the idealised data set, so its about 1000 sequences and the tree will guide you into major groups - which is probably what is wanted, i.e. to place a 16S from a local salmon farm into the broader clade diversity of the species.

---

The answer is for this species is there are around a thousand sequences on Genbank for the species of interest which are split into 6 identical genotypic groups. The demonstration of the species is represented by the scale bar, i.e. it is very little genetic diversity.

*Cycloclausticus sp* is the outgroup and each member of the in-group is *Piscirickettsia salmonis*

---

Just to address the comments

* "Cycloclasticus"? The question is not specific to /db\_xref="taxon:1238"

It is specific there's a fundamental principle of trees. A tree needs a root, all the ingroups are all *Piscirickettsia salmonis*. Without a root, i.e. something that is not *Piscirickettsia salmonis* there is no direction in the tree. *Cycloclasticus* is a good or reasonable root.

* how do i access this interface?

Its complicated to explain and its a different question.
> -1 votes

---
Tags: genbank, sequence, 16rrna, ncbi
---


title: How do I rename fasta headers?

text:
# Question
Title: How do I rename fasta headers?
I am trying to use the seqkit replace command to replace chromosome names in the format `chr_I, chr_II, ...` to `I, II, ...`. I am using the following command:

```
seqkit replace -p "(.)" --replacement "{kv}" --kv-file keyValues.txt mySequence.fasta

```

My `keyValues.txt` file contains the following:

```
chr_I   I
chr_II  II
...

```

The two columns are separated by a tabulation.

I get the following output:

```
[INFO] read key-value file: keyValues.txt
[INFO] 6 pairs of key-value loaded
[ERRO] pattern "(.)" matches multiple targets in "chr_I", this will cause chaos

```

I don't know what I am doing wrong. I chose the `(.)` pattern to match the whole header but it seems to be wrong. Any help would be appreciated.

### Update

All the headers are shown below:

```
grep '>' mySequence.fasta
>chr_I
>chr_II
>chr_III
>chr_IV
>chr_V
>chr_X

```

# Answer
> 10 votes
Why not use `sed` instead?

```
sed -e 's/chr_I/I/' -e 's/chr_V/V/' -e 's/chr_X/X/' mySequence.fasta > mySeq.fasta

```

Or even simpler:

```
sed 's/chr_//' mySequence.fasta > mySeq.fasta

```

# Answer
> 2 votes
I like using `awk` when there's any kind of lookup involved: something like

```
awk 'FNR==NR {
    hash[">"$1]=$2;
    next;
}

if ($1 ~ /^>/) {
    print ">"hash[$1];
} else {
    print $1;
}' keyValues.txt mySequence.fasta

```

Basically the `FNR==NR` check tells awk to work only on the first file (useful primer here), `keyValues.txt` and create an association for each key (eg. `chr_I`) with its value (eg. `I`). The rest of the code after the `next` works only on `mySequence.fasta`, printing out the lookup value only if the line is a fasta header, as checked by the `$1 ~ /^>/` condition.

# Answer
> 1 votes
Here is an alternative to sed: `cat test | tr -d "chr_"`

The pipe is to output the file, while the `tr` command is to translate and change the file

# Answer
> 1 votes
Your issue is that the pattern `.` matches a single character. To match the whole header, use `.*`:

```
seqkit replace -p "(.*)" --replacement "{kv}" --kv-file keyValues.txt mySequence.fasta

```

Alternatively as pointed out by others you may not really need `keyValues.txt` here as you're just removing the `chr_` prefix:

```
seqkit replace -p "chr_(.*)" --replacement '$1' mySequence.fasta

```

(Don't forget the single quotes for `$1` else the shell will expand it)

---
Tags: fasta
---


# Question
Title: How do I save selected sequences in seqkit to a file?
Consider:

> seqkit grep -nrp ";s\__nosocomialis;" ribogrove\_11.217\_sequences.fasta

Using this command lets me select all sequences related to the organism "*A. nosocomialis*" in the RiboGrove database.

How do I save these selected sequences to a new ".fasta" file?

# Answer
Based on the docs, seqkit grep has a `-o` "output file" option, so if you get the correct output with:

```
seqkit grep -nrp ";s__nosocomialis;" ribogrove_11.217_sequences.fasta

```

Then this command will (in theory) give you your desired outcome:

```
seqkit grep -nrp ";s__nosocomialis;" ribogrove_11.217_sequences.fasta -o new.fasta

```
> 4 votes

# Answer
For any program that prints output to the terminal (standard output), you can always redirect to a file using `>` or `>>`. Either of the two will create the output file if it doesn't exist, but `>` will overwrite any contents of an already existing file, while `>>` will append to it without overwriting. So you could just do:

```
seqkit grep -nrp ";s__nosocomialis;" ribogrove_11.217_sequences.fasta > out.fa

```
> 1 votes

---
Tags: seqkit
---


# Question
Title: Merging multiple seurat datalist together
I have a R script that subsets out the tumor cells from 6 different individual seurat objects.

It gets a list of file paths for all RDS files containing the expression data for multiple datasets in a directory. It creates two empty lists, datasets and dataset\_index, to hold the subsets of tumor cells from each dataset and their corresponding indices, respectively. It loops through each file path, reads the RDS file, and subsets the tumor cells from the expression data. If the dataset does not contain any tumor cells, it prints a message indicating that fact. The subsets of tumor cells are appended to the datasets list and their corresponding indices are appended to the dataset\_index list.

```
scfiles = list.files('/LungDataSets/',pattern = '.*SCTransform_Symphony.*.RDS',recursive=TRUE)

scfiles = paste('/LungDataSets/',scfiles,sep='')

datasets = list()
dataset_index = list()

for (i in 1:length(scfiles)) {
    scData = readRDS(scfiles[i])

    if (is.element('Tumor', unique(scData$cell_ann)))  {
        tumorData = subset(scData,subset = cell_ann=='Tumor')
        datasets = append(datasets,tumorData)
        dataset_index = append(dataset_index,toString(i))
    }
    else {
        print(c(scfiles[i], 'does not contain Tumor cells'))
    }

    rm(scData)

}

```

Now I want to merge all those tumor cells from the 6 datasets into one seurat object.

I've tried `merge` and `rbind` but they don't work.

```
megaObject <- merge(datasets)
Error in as.data.frame.default(x[[i]], optional = TRUE) :
  cannot coerce class ‘structure("Seurat", package = "SeuratObject")’ to a data.frame

```

# Answer
The first parameter of `merge` should be a Seurat object, the second (`y`) can be one Seurat object or a list of several. Try:

```
merge(x = datasets[[1]], y = datasets[-1])

```

See the merge vignette for more details. You may want to use the `add.cell.ids` option to be able to tell which dataset each cell originated from.
> 3 votes

# Answer
You should be able to use something like this:

```
Reduce(
    f = function(x, y) {merge(x, y, merge.data = FALSE)},
    x = datasets # list of Seurat objects
)

```

This will create a new Seurat object based on the multiple seurat objects in your list. It might be good idea to store the "sample" information within the metadata slots of individual objects. Probably a better approach is to have cell barcodes of each sample prefixed with the sample id.
> 2 votes

# Answer
I can't reproduce things quickly without a minimum reproducible example (e.g. with sample data from the package). However, the error is telling you that R cannot convert the object to a dataframe when trying to use `merge()` (at least I suspect that this error occurred when you ran `merge()` since you posted it right under it). I suggest you convert *everything you need* to a dataframe (there are quite some questions, e.g., on biostars regarding this conversion). *Everything you need* because a SeuratObject cannot simply be cast to a dataframe, only parts of it, e.g., `data_transcripts_seurat@assays$SCT@counts %>% as.matrix() %>% t() %>% as.data.frame()` (from this biostars question). However, you might need some other information. Therefore, I suggest you explore the SeuratObject using RStudio and check what information you want to draw from it. After the (interesting) parts of the SeuratObject have been cast to dataframes, e.g., in the for-loop, you can go ahead with merging.
> 0 votes

---
Tags: r, scrnaseq, seurat, merge
---


# Question
Title: GATK SplitIntervals: A USER ERROR has occurred: but no positional argument is defined for this tool
I am trying to create a GATK g CNV pipeline and I need to make shards before running the `GermlineCNVCaller` therefore, I believe I need to use `SplitIntervals`.

I have the following predefined variables:

* `singularity_image` is the path to the GATK sif file
* `bam_file` is a bam file from a defined bam directory
* `reference` is a path to a .fa file
* `targets_interval_list` is a interval list created by Picard (ending with .interval\_list)
* `NUM_BINS` is a manually preset value from 100 to 200.

```
  singularity exec $singularity_image /gatk/gatk --java-options -Xms8000m \
  SplitIntervals \
  -I ${bam_file} \
  -R ${reference} \
  -L ${targets_interval_list} \
  --scatter-count ${NUM_BINS} \
  -O ${output_dir}/Interval_list_${bam_filename}_.shard_\%${NUM_BINS}\%_\%{shard_index}\%.interval_list

```

I have tried change the specifications under `SplitIntervals` but for whatever reason I keep getting the following error

```
***********************************************************************

A USER ERROR has occurred: but no positional argument is defined for this tool.

***********************************************************************

```

Why is this error occurring?

**The documentation I am using is here:** https://gatk.broadinstitute.org/hc/en-us/articles/360036899592-SplitIntervals

# Answer
> 1 votes
This solution where the `-O` is just set to the output directory seems the most easy way to handle the issue:

Keep the following variables the same:

* `singularity_image` is the path to the GATK sif file
* `bam_file` is a bam file from a defined bam directory
* `reference` is a path to a .fa file
* `targets_interval_list` is a interval list created by Picard (ending with .interval\_list)
* `NUM_BINS` is a manually preset value from 100 to 200.

The code with the modified `-O`:

```
  # run SplitIntervals with shard name
  singularity exec $singularity_image /gatk/gatk --java-options -Xms8000m \
  SplitIntervals \
  -R ${reference} \
  -L ${targets_interval_list} \
  --scatter-count ${NUM_BINS} \
  -O ${output_dir}

```

While it doesn't name the shards in a very elegant way it seems to work. The shards are made from the interval list and therefore, should not be named after the .bams that will be processed in other steps that are used (for in this case a GATK g CNV pipeline).

---
Tags: gatk, singularity
---


# Question
Title: What is the limit of partition incompleteness in Bayesian MCMC?
I am seeking advice on whether or not to incorporate incomplete data into my **BEAST** analysis.

I have HIV data consisting of two partitions of the *pol* gene: one that is 1.2 *kb* in length and present in all isolates, & another that is 2 *kb* in length but only present in one-third of the isolates. I am wondering if I should include the second partition in my analysis, despite its incomplete nature.

I have found some publications discussing the use of incomplete data in phylogenetic analyses, such as Yuchi Zheng & John J. Wiens (2015), Rasmus Hovmöller et al. (2013), John J. Wiens & Jonathan Tiu (2012). However, I remain uncertain about the best course of action.

If you have experience working with incomplete data in BEAST analyses or any insights on this matter, I would greatly appreciate your advice.

Thank you in advance for your help.

Best regards,

---

## Update

As seen from the answer & comments, there has been a misinterpretation of my partitions' structure & their incompleteness. So hereinbelow is a short description of my data, the subject topic, & a diagram.

There are 1009 isolates of the *pol* gene of HIV from the Republic of Armenia. They were sequenced from the start of the gene by two means: Sanger & an NGS; which caused two groups of length: 1.2 kb & \> 3 kb, respectively.

This way the situation may be demonstrated via the next artificial miniature diagram:

Based on the picture like above, I have thought it rational to divide the data into two staple partitions: the first third of *pol* & its last part. The first partition has all sequences completed, and the second one has mostly full-length complete sequences for a third of isolates and just `NNN…NNN`s for the other two-thirds.

`N`s are automatically added at the alignment step. Sanger sequences are just short and end from 1200 to 1500. I currently do not have access to the source data. They were provided to me already built and aligned at the reference.

In summary, the question is, **"Should only the complete first partition be used in subsequent BEAST analysis or should the incomplete second one be included as well?"**

# Answer
> 1 votes
The original answer was "criticised" for being obtuse.

The basic very simple answer is:

1. If the data you present in the diagram (alignment) is left EXACTLY in the format used to make the tree - it honestly does not matter whether you have 1 partition or 2 partitions. It makes no practical difference. The posterior likelihood and stuff will change a bit - basically the result is the same. The reason is given in the very long explanation below. I admit 2-partitions looks cooler and there's reason not to do it, computationally it will be a bit slower. **Partitions** I am happy to explain why a partition is used - but thats different answer or question.
2. What must not be done with HIV is changing every '-' to '?' or changing every '-' to 'N'. Whether this is 1-partition or 2-partitions. With higher eukaryotes that is fine and that will work very nicely, but HIV no. So Weins et al (2012) can do this (almost certainly have done), but you can't.

I know it seems a very small change replacing every '-' with a '?' or 'N', but what it represents seismic. I have attempted to explain why it is so important below the line, but it's about indels (insertions/deletions). I recognise my explanation may not be perfect.

**Confusion** The confusing bit is you mention 'NNNNNs' and thats the worry. If NNNNs are changed to '-' thats fine, if they are left as 'NNNNs' that raises the complicated spectre of 'indels' and will get very complicated. It probably can be done but its really complicated.

To make the point clear if I perform a `trimal` analysis on my alignment

```
trimal -in infile -out outfile -nogap -fasta

```

It will remove all sites with '-' in them, i.e. it will remove the entire site, not just the '-' character. I think it's really cool. Anyway, the tree from a `trimal` alignment will be basically the same as '1-partition' or '2-partitions' with all the '-' retained, even though the alignments look very different. That sounds weird but it is a consequence of how Bayes (i.e. Beast) or ML calculate phylogenetic information, which I tried to explain below. HOWEVER, if there is an alignment with 'NNNNs' instead of '-' that will produce a very different tree.

So for simplicity the answer is that 2-partitions is fine for '-', but not if those '-' are represented as 'NNNNN'. Thus for '-' alignments it looks cool, but makes little difference. 2-partitions can make a HUGE difference with different data particularly with loads and loads of NNNNNs instead of '-'.

---

**Summary** HIV is a special case because *env* in particular - and possibly all the HIV genes - are known to contain extensive indels. Indels in *env* are so prevalent in HIV they can be approximated via distribution (did Ziheng Yang do that?). So the short answer, Weins *et al* (2012) can do it because they're dealing with higher eukaryotes and genes with negligible indels (in addition they could remove known indel sites), in HIV you would need independent evidence the missing data site was not going to be an indel (i.e. that particular site has never had a known indel in any HIV data set).

The crux is there is a difference between missing data and absent data. Absent data *can be* a no-go (below; context depending - some times it makes no difference \[immediately below\]), it can be detrimental to the phylogenetic signal. Missing data is very interesting.

Basically, if you are going to leave the missing 2kb gene as '-' ... thats perfectly fine. It will not add any phylogenetic information - but it will not hurt the calculation. If 'N' or '?' are used instead to represent a missing 2Kb locus for HIV instead of '-' thats different.

---

Weins *et al* (2012) is a very interesting study.

There are four issues

1. Sample skews
2. Absent data which causes blockwise deletions in Bayes and ML
3. Missing data which *might* recover phylogenetic information at a given site if correctly annotated. 4.Quality of data

The traditional view is that it is better to remove the site/ or sample, if the quality of data is in doubt or for example an entire gene is missing. The traditional view is it is better to be certain of the data regardless of possible risk of a sample skew. However, *sample skews* are the Achilles heal of population genetics and a lesser extent to phylogenetics. Thats what Weins *et al* (2012) is concluding, it's really cool. I like it.

We all know that an analyst can't deliberately cause a sample skew (i.e. you can't omit samples without good reason) - that goes with out saying. Right?

The context of the question is "if a sample has *pol* and its other gene wasn't sequenced - do you fill in the absent gene with "Ns" or "?" to trap the phylogenetic information in *pol* and the samples which have data in the gene/locus/partition with limited coverage. **It's a really good question.** My suspicion for HIV is I'm afraid you probably can't, but that might just be *env*.

**Absent data** Bayesian and ML have no '-' character in their matrix. If an indel is present it will trigger that site to be excluded from carrying any phylogenetic information specifically for Bayes and ML. Thus

```
                Site1
Sample 1         A
Sample 2         A
Sample 3         A
....
Sample 9999      -
Sample 10000     A

```

Site 1 will not carry any phylogenetic information for Bayes or ML because the single '-' is not represented in the mutation matrix. The entire site for 9999 nucleotides is excluded. A pairwise analysis would be needed, e.g. distance methods.

**Missing data** However "Sample 9999, Site 1" could be represented as 'N' instead of a '-' because that means any nucleotide in the matrix. *That works!* Site 1 is now phylogenetically informative under Bayes and ML. Hooray! However, are we *certain* Sample 9999, Site 1 is a nucleotide rather than an indel? If indels are rewritten as nucleotides, i.e. an absent nucleotide is converted to a nucleotide of unknown type ('AGCT'), i.e. missing data - non-nucleotide to nucleotide - thats cheating.

This creates issues **particularly** for HIV *env* where there are loads of indels. If an indel was replaced with an 'N' - the analysis might work nicely, but there is real risk it is not being true to the data.

---

**Possibly explanation** Sample skews in phylogeny are really interesting. What can happen is long internal branches can cause changes to the branching order. However, if those internal branches are shorten by more rigorous taxon sampling then these shortened internal branches can recover the true true. I think it was Jo Felsenstein that discovered this (the guy who invented maximum likelihood as a method for phylogenetics in place of parsimony).

One explanation of Weins *et al* (2012) is that the 'missing data' samples are shortening the distance between internal branches - by providing a more rigorous sample - and that is resulting in recovery of the correct tree. So it is probably a result inline with Jo Felsensteins original observation.

So if ***pol*** never has indels - ever - you can do it. This would also help population dynamics and all shorts of calculations. If there's a risk that you are replacing indels with nucleotides it is very much, much better to avoid it.

**Conclusion** Everyone in phylogentics knows the issues of sample skews - a possible sample skew is permitted if the data is genuinely never been produced - it's not in existence (obviously). However, it is not permitted if the data exists but the analyst hasn't included it - for whatever reason. Its also preferred if the quality is in question - its better to err on the side of caution.

If you fill in an absent gene with "?" or "N" to obtain phylogenetic information from genes that are present for that sample ... that is permitted providing it doesn't change the possible nature of the underlying data. *The rationale for doing this is the risk of sample skews.*

---

So would I adopt a "non-traditional" approach? As a general rule for viruses no and this is a particular special case regarding indels.

If this was a higher eukaryote and like - that critical sample and *that* gene - didn't give sequence. This has happened. So would I dump the sample (or in this case the locus) - no I'd caveat the tree and explain and empirically check the assumptions.

If I was in a situation where 2/3 of the data were missing, for a higher eukaryote (DNA with great proof-reading), personally I would do a full scale investigation, i.e. phylogenies with and without the gene/locus/partition in question. That is pretty much was Weins et al (2012) have done.

---
Tags: phylogenetics, phylogeny, virus, beast, hiv
---


# Question
Title: Why is there antisense sequence in RNAseq data
I'm looking at RNAseq data from CCLE. The data is paired-end.

Take the cell line Hs578T and the gene *HRAS* as an example.

The cell line carries a G12D mutation (c.35G\>A), so the change in cds is:

```
ggc ggtgtgggca agagtgcgct g - Wildtype CDS
gAc ggtgtgggca agagtgcgct g - Mutant CDS
 ^

```

My question is, when I `grep` the mutant CDS `gAcggtgtgggcaagagtgcgctg`, I do not get a match in my .bam file. But when I `grep` the reverse complement (i.e., anti-sense sequence), I get matches coming from both mates.

My question is, why are there reads in the .bam file that correspond to the anti-sense strand of the DNA? Shouldn't all the reads of mRNA be the CDS?

---

I got the answer - The punch line is that the sequence in .bam files is ALWAYS '+' strand of the reference, no matter what mate the read came from, or which strand the gene sits on.

# Answer
> 3 votes
Yes, the mRNA results should include CDSs (not only, don't forget you also have UTR sequences which are non-coding so are not part of the CDS), but there is no reason why they should all be from genes on the + strand of the genome. Remember that the distinction between "sense/+" and "antisense/-" is essentially a human invention; it doesn't imply that the cell considers one strand to be in any way superior than the other, and genes are pretty randomly distributed across both strands.

This means that when you align to a human reference genome, you will always have some sequences aligning to the + strand, if the RNAs were transcribed from genes that happen to be on the sense strand, and others will align to the - strand, if they were transcribed from genes found in the antisense strand.

This is completely normal and expected, there is nothing wrong with your data. For example, human tP53 is on the antisense strand (note how the arrows are all pointing to the left):

---
Tags: rna-seq, strandedness
---


# Question
Title: How to analyze hashtagged data with scanpy?
I used hashsolo to demultiplex my data and I am a bit confused as to how to continue the analysis. Following hashsolo my samples look like below:

```
A1
AnnData object with n_obs × n_vars = 8684 × 18197
    obs: 'most_likely_hypothesis', 'Classification', 'n_genes'
    var: 'gene_ids', 'feature_types', 'genome', 'n_counts', 'n_cells'
    layers: 'counts'

A2
AnnData object with n_obs × n_vars = 9901 × 20862
    obs: 'most_likely_hypothesis', 'Classification'
    var: 'gene_ids', 'feature_types', 'genome', 'n_counts'
    layers: 'counts'

```

I have a few questions regarding the analysis. My initial plan was to merge the 2 samples. However, when I do that and I print the Classification the merged object seems to add the classification counts because the same hashtag names were used.

For example hashtag 1 has 949 in sample A1 and 1286 in sample A2. When I do print`(adata.obs['Classification'].value_counts())` hashtag1 has 2235 counts. Will that be a problem with the analysis? I merged the object like this:

```
# add some metadata
A1.obs['type']="C"
A2.obs['type']="T"

# merge into one object.
adata = A1.concatenate(A2)

adata
AnnData object with n_obs × n_vars = 18585 × 18064
    obs: 'most_likely_hypothesis', 'Classification', 'n_genes', 'type', 'batch'
    var: 'gene_ids', 'feature_types', 'genome', 'n_counts-0', 'n_cells-0', 'n_counts-1'
    layers: 'counts'

print(adata.obs['Classification'].value_counts())
Hashtag_2_Antibody    5134
Hashtag_3_Antibody    3150
Hashtag_4_Antibody    2000
Hashtag_1_Antibody    2235
Doublet                                         1598
Negative                                           9
Name: Classification, dtype: int64

```

My hashtags also show up wen I do `sc.pl.highest_expr_genes(adata, n_top=20, )`. Is that a problem?

I would also like to plot a UMAP following leiden for each hashtagged sample but I am a bit confused on how to do that. Normally, I would add information such as sample or type on metadata and color by sample for example. But if I try to color by classification it will merge the hashtags from both samples.

Thank you

# Answer
> 0 votes
You could add something like `A1.obs["Classification_obj"] = [c+"_A1" for c in adata.obs["Classification"]]` but the classification should also be clear in the merged adata object from the "Classification" and "batch" labels together

More importantly, after running hashsolo, you must delete your hashtag feature\_types from adata.var, something like `A1 = A1[:, [g for g in A1.var.index if g not in list_of_Ab_hashtags]]` . Once you have the cell type classifications and are happy with them (can be worth doing further QC first), you can throw away that Ab count data for the downstream standard GEX analysis

---
Tags: rna, umap, hashing
---


title: mapping nucleotide sequences

text:
# Question
Title: mapping nucleotide sequences
I am trying to do mapping with multiple sequences against a reference genome. I tried the following code: The outputfile is .fastq.sam. I need only .sam. How do I get it?

```
for i in test/*.fastq;
do
 ./minimap2 -ax map-ont reference.fna $i > $i.sam;
 done

```

# Answer
The command is fine, it is just a suffix. If you want to trim the suffix then use `${i%.fastq}.sam` for the output, or even better convert the output from minimap2 to bam directly as this saves space:

```
./minimap2 -ax map-ont reference.fna "${i}" | samtools view -o "${i%.fastq}".bam

```
> 5 votes

# Answer
Honestly, I wouldn't do it with a one-liner like this. Make a proper bash script, so you can spend a couple of lines to get the naming just right. Then you can use the bash script over and over again. Easier to document too if the command line used is in a file, not just in terminal history.

Also, yes, piping to convert to bam is almost certainly what you want to do. But piping into a sort commmand is better. And I think samtools sort will output a bam by itself if you ask it to.
> 1 votes

# Answer
I would recommend

```
ls test/*.fastq|sed s,.fastq$,,|xargs -i echo ./minimap2 -ax map-ont reference.fna {}.fastq \> {}.sam|sh

```

Replace trailing "sh" with GNU parallel for concurrent jobs on the same node; replace with asub for submitting array jobs to Slurm/LSF/SGE.

Once you get used to this, you can do a variety of tasks quickly, like renaming files in batch, gzipping all files in parallel, format conversion in batch, aligning 100 files on a cluster, etc.
> 1 votes

---
Tags: read-mapping
---


# Question
Title: Gene rank scores from LINCS database: Evangelista et al. 2022
The website https://maayanlab.cloud/sigcom-lincs/#/SignatureSearch/UpDown can be used to "Identify reversers and mimickers from over 1 million signatures by entering up and down gene sets .. " I am looking to get these "reversers" and "mimickers", aka chemical perturbagens.

So on this website, I paste my list of a selection of Up-regulated genes, in their respective box, and the same for my list of down-regulated genes. These genes are from an RNA sequencing experiment.

After the website has done its thing, I expand the search for chemical perturbagens (LINCS L1000 Chemical Perturbations (2021)). One of the data tables I can download contain something called a Gene Rank Matrix. This table consist of one row for each gene and columns of the chemical perturbagens found. The values are integers.

A small subsection of the data table would look like this:

My questions are what these integers (12, 144, 411, 641) represent. Are they how much the perturbagens (In this case REP.B025\_YAPC\_24H\_K01\_resiquimod\_2.22uM) will affect the genes?

The paper describing this website searches (https://academic.oup.com/nar/article/50/W1/W697/6582159) mentions "The expression tables provide the computed CD differential gene expression coefficients for each signature, while the rank matrices provide the integer gene ranks for each signature as determined by the coefficients"

and

"Matrices that represent the gene rankings of the L1000 signatures are stored as HashMaps by the data API in random access memory (RAM). HashMaps are data structures designed to efficiently search for key-value pairs, and this implementation drastically cuts the fetch times of the gene rankings by UUIDs, thus improving the speed of the search"

Does this mean that the scores found wont have any meaning in regards to how well the perturbagens found will affect the genes? Or can some meaning be derived from this?

Thanks in advance!

# Answer
**Summary answer** Although I'm unhappy about how the paper is written and possibly the methodology, the following sentence is likely (almost certainly) to be the answer.

The coefficient (aka weight) is how important that gene is for the classification of that signature. Thats it! :-)

Essentially, there are loads of genes involved in identifying which genes represents which signature. Not all those genes are equally informative. The positive or negative skews in that calculation, called a "classification", are identified through their weights (aka coefficients). Thats simply the answer. The bigger the coefficient the more important that gene is in the classification, i.e. its associated with the signature.

There is no need to be more complicated than that. If you are interested in a critique of that paper, which I had to do to get to the answer and thereby understand the analytics of that paper, this is presented below.

---

---

**Critique of Evangelista *et al*. 2022: Nucleic Acids 50:W697–W709**

**Note 1** I don't have time to read the paper in absolute detail, particularly the appendices, where some of the information needed may have be placed.

**Note 2** In the answer below, a "expression coefficient" is what I've called an "expression weight" its the same thing.

**Weightings and analytical concerns** I do understand weightings in deep learning, particularly machine learning. I read the paper and the ANN (artificial neural network) was very basic. I don't believe it had multiple layers, because it would have said so and it was a **very odd** description for a deep learning algorithm. It would also be odd to have a single layer in deep learning, *very odd*. The loss function is also not described at all which is *extremely important* in deep learning (it will not work otherwise). Essentially, what the authors call "deep learning" appears to a machine learning neural network using an Adams gradient descent. Adams gradient descent is a standard approach in deep-learning admittedly. Thereafter in the paper the details were very sketchy, maybe they are in the appendix (should be checked). Any deep-learning, machine-learning requires a demonstration of the robustness of the training and this appears absent. The two basic stats are accuracy and AUC-ROC, which were definitely absent. Universal stats are used to quickly understand how good the algorithm is in comparison to other deep-,machine- learning algorithms for other data sets.

**Analytical overview** Generally I would assume the authors have performed a supervised learning algorithm, but this may not be correct (below). If it is supervised learning they have performed they would have identified their 'signatures' against a very large data set, then trained the different types of signature against that data and tested it. The learning could be terrible - but it's just not stated at all. I'm surprised because Nucleic Acids is at the pinnacle of algorithm publications.

**Summary.** The authors performed deep-learning (which looked like a neural network machine learning algorithm). The description is at best poor, but further information could be present in the appendix. Therefore its not absolutely clear what a "weight" would represent because there are two possibilities.

**Weights possibility: unsupervised learning** The one thing that might have happened is they've performed an unsupervised learning algorithm and the weights are a sort of "maximum variance". This very trendy. That would explain the lack of accuracy and AUC-ROC type stuff. **IF** this is true and I don't believe it's explained - I could be wrong - then the weights represent **maximum information**. What does that mean? Good question, its context specific, its the gene which is most important in sectioning up the data - e.g. housekeeping gene which is required 'all over the place', not just part of the cell cycle is an example. Basically, in unsupervised weighted output high ranking weights are a mix of really boring and really brilliant stuff. So you just go through it with biological understanding, is this gene boring yeah/neigh, is this one - wait wow! unexpected ... That is actually fashionable analytics :-)

**Weights possibility: supervised learning** Is already answered at the beginning. This is almost certainly the correct answer - it would be the orthodox approach - and we just accept the robustness checks for the algorithm and possibly the algorithm itself were simply insufficient.

**Why isn't the paper clear (at best)?** What I assume has happened is the team have brilliant frontend and genomics team backend. The analyst in the middle is unlikely to be a specialist and thats why it's poorly described, or possibly performed. I assume if the reviewers were of similar composition then that's how its got through Nucleic Acids. Thus, I would personally treat the output with caution because deep learning (machine-learning) will be the essential bottleneck between their backend development and your data. Maybe "bottleneck" is not quite an accurate description, because its speeding up the classification, so maybe a "gate-keeper" is a better description.
> 1 votes

---
Tags: rna-seq, public-databases, machine-learning, drugs, deep-learning
---