# Question
Title: Can not launch bcftools using python's subprocess module, as it only accepts first command of commands list
I am trying to remove samples from a chromosome vcf file. I wrote a function that takes chromosome number and a list of samples to remove. When I try to run bcftools using subprocess module it only runs bcftools, as if I was running `subprocess.run([bcftools], shell=True, env=my_env)` instead of:

```
subprocess.run(['bcftools', 'view', in_vcf, '--samples', drop_samples, '-Oz', '>', out_vcf], shell=True, env=my_env)

```

The function:

```
def remove_individuals(chr_samples):
    my_env = os.environ.copy()
    my_env['PATH'] = f'/home/aydar/bin/bcftools_1.11/bin/:{my_env["PATH"]}'
    def remove_indiv(chr, indiv):
        in_vcf = f'data.filtered-updated-{chr}.vcf.gz'
        out_vcf = f'data.filtered-updated-{chr}-faulty-samp-rem.vcf.gz'
        drop_samples = "\'^"+",".join(indiv)+"\'"
        proc = subprocess.run(
            ['bcftools', 'view', in_vcf,
             '--samples', drop_samples, '-Oz', '>', out_vcf],
            shell=True,
            env=my_env
        )
        print(proc.args)
    chr_samples.apply(lambda x: remove_indiv(x.index[0], list(x.values[0])))

```

function input:

```
13,[2500000075_2500000075,2000000171_2000000171,1800000179_1800000179,330134048_330134048]

```

When I run the same command in the terminal everything works.

The command I test with in terminal:

```
bcftools view data.filtered-updated-13.vcf.gz --samples '^2500000075_2500000075,2000000171_2000000171,1800000179_1800000179,330134048_330134048' -Oz > out_test_vcf.gz

```

# Answer
> 3 votes
Never mind, I figured it out. As it turned out I had to write to file like so and extra quotes `("\'")` around samples were not needed. Also shell=True does not work well with lists, only strings:

```
def remove_individuals(chr_samples):
    my_env = os.environ.copy()
    my_env['PATH'] = f'/home/aydar/bin/bcftools_1.11/bin/:{my_env["PATH"]}'
    def remove_indiv(chr, indiv):
        in_vcf = f'data.filtered-updated-{chr}.vcf.gz'
        out_vcf = f'data.filtered-updated-{chr}-faulty-samp-rem.vcf.gz'
        drop_samples = "^"+",".join(indiv)
        out = subprocess.run(
            ['bcftools', 'view', in_vcf,
             '--samples', drop_samples, '-Oz'],
            env=my_env,
            capture_output=True
        )
        with open(out_vcf, 'wb') as f:
            f.write(out.stdout)
    chr_samples.apply(lambda x: remove_indiv(x.index[0], list(x.values[0])))

```

Remember to use `'wb'` if you like me and use `Oz` output option, because it returns bgzipped file.

---
Tags: python, bcftools
---


# Question
Title: Variation in 16S rRNA between assemblers - how do I know which is real?
I have a low-diversity metagenome (~11 bins \> 80% completion). Out of the bins, 3 are of interest to me. None of the lower-completion bins that can be identified are from the group of interest. So let's call the bins of interest 1, 2 and 3.

The raw reads had been assembled separately with MEGAHIT and SPAdes (meta). Out of them, only the SPAdes assembly was binned / refined. There are 3 16S rRNA in the MEGAHIT assembly (x, y, z) and 3 16S rRNA in the SPAdes assembly (a, b, c).

Bin 1 is associated with sequence 'a', which is identical to 'x'.

Bin 2 is linked with SPAdes sequence 'b', which has no MEGAHIT equivalent. Bin 3 is not associated with any 16S rRNA (even with MarkerMAG).

Sequences 'c' and 'z' are very similar. Sequence 'y' is an outlier.

There are small hints from phylogenetics that bin 3 could be associated with sequences c/z. But then, could sequence 'y' just be a misassembly? Or is it that somehow my final bins are missing an entire genome of interest? If the 16S rRNA for 'y' is real, I would expect a genome divergent enough from the others that it wouldn't get binned together with them. This is what I am trying to understand. Would you happen to have any advice / suggestions on what I could to to test this?

Initial approach: assessing scaffold quality with the BBMap stats.sh . It's a bit difficult to tell, since SPAdes has about 10x more contigs. Worse N50/N90 but more long contigs and a better 'general assembly score'.

Additional thoughts: The new bins seem too divergent from known genomes to use a reference for misassembly estimation.

Thank you for your time!

# Answer
This is easy @Laura - I'm a bit surprised by the results Spades meta didn't have a good reputation, but the criteria assessed was not what you are looking for.

Anyway take all outputs and place them in large multi-species alignment,

* e.g. for each output perform a Blast for maximum genetic diversity for each sequence and sample across the large Blast output. Take note if there are multiple partial hits against diverse bacteria (online Blast has a diagrammatical summary).
* You have a good understanding of your taxonomic group and intersperse that with other stuff, i.e. fill in the diversity.
* Assemble it, e.g. muscle5 superfast option.
* Perform a split decomposition phylogeny or equivalent which identifies phylogenetic networks.

SplitsTree is the historic package for this and it turns out that split decomposition was not the best, I think "neighornet" was superior. What this is does is identify mosaic gene signal in your data. It provides a visual representation. Thus if your data is a artifactual mix of two or more taxa this will cause a network effect between the two parents and your newly mis-assembled sequence will sit between the two in a two "diamond same". This is the antonym of a bifurcating phylogeny.

The final step is the assembly that produces sequences that minimises phylogenetic networks is the true sequence. If you've not encountered split-decomposition/networks before a verbal explanation might seem confusing but when you run through a few examples and try it on your data you'll get the grasp.
> 2 votes

---
Tags: phylogeny, software-recommendation, metagenome, 16rrna
---


title: Singularity, AnnotSV and HPO

text:
# Question
Title: Singularity, AnnotSV and HPO
I am trying to use AnnotSV to annotate some Manta outputs in the following code with the HPO term changed:

Where:

* input\_dir is the input directory
* singularity\_container is the path to the AnnotSV singularity container

```
# Loop through the files in the input directory
for file in "$input_dir"/*.vcf; do
    # Get the base name of the file without extension
    filename=$(basename "$file" .vcf)

    # Create the output subdirectory within the input directory
    output_dir="$input_dir/AnnotSV_Batch_Run"
mkdir -p "$output_dir"

    # Set the output file path with "_AnnotSV" appended
    output_file="$output_dir/${filename}_AnnotSV_"

    # Run AnnotSV via Singularity with the base name as output prefix
    singularity exec "$singularity_container" AnnotSV -SVinputFile "$file" -genomeBuild GRCh38 --hpo "HP:0000001"  -outputFile "$output_file"
done

```

The HPO term I am using is in the jax.org site but when I run it - I get the most peculiar error  `option not known.` but the HPO term is specified in the AnnotSV manual on page 49 (according to the number on the bottom). As for trying I tried with `--HPOterm` and `hpoterm` I tried with and without the `:` and with and without quotes as well. I am not sure what the issue is but I cannot seem to find the solution. Has anyone perchance used this software before with HPO terms as it seems not to be accepting the input I give it?

# Answer
> 1 votes
It seems that currently AnnotSV does not support HPO terms usage an Singularity.

---
Tags: singularity
---


# Question
Title: Snakemake: Re-read the samples.tsv file midway through the workflow
All, I have a workflow that in some rules requires the sample's read-lengths. For this I always ran a script that would parse all fastq files and would add the read lengths to the sample table. Rules could then take the read lengths for each sample from there during execution.

However, I recently made a workflow/pipeline that downloads files from SRA, so I need to compute the read lengths at a later stage during workflow execution (after prefetch/fasterq-dump).

I added a step to update the sample table, however, I think the sample table is only read at the very beginning of the workflow.

What would be the nicest way, the Snakemake preferred way, to re-read the samples.tsv file before a rule is executed?

In general I like storing meta data in the sample table so, I'll probably use this mechanism more often in the future, if I find a nice methods to re-parse the samples table.

Right now I do this in rules:

```
rule deletion_profile:
    input:
        rseqc_input_bam,
        rules.add_max_read_length.output.done_file
    output:
        os.path.join(rseqc_dir, '{sample}.deletion_profile.txt')
    conda:
        "../envs/rseqc.yaml"
    params:
        read_length = lambda wildcards: samples.loc[wildcards.sample, 'read_length']
    shell:
        '''
        deletion_profile.py -i {input} -l {params.read_length} -o {rseqc_dir}/{wildcards.sample}

```

Thanx.

# Answer
> 3 votes
You're right, it doesn't re-run things again for each job, so your original definition of `samples` is still the one it's working with.

It seems like it goes against the grain of Snakemake to overwrite a file in-place that acts as input to a rule, so I think the most Snakemake-y way you could do it would be to treat your samples.tsv as part of the input/output graph like everything else. You could have a rule with samples.tsv as input and samples.somethingelse.tsv as output where that latter one has read lengths and whatever else gets automatically added. (Alternatively, you could just have your updated samples.tsv read within the rule itself rather than relying on a shared copy loaded earlier, which is kind of halfway to the above idea.)

Either way, maybe something like:

```
rule deletion_profile:
    input:
        bam=rseqc_input_bam,
        samples="samples.extended.tsv"
    output:
        os.path.join(rseqc_dir, '{sample}.deletion_profile.txt')
    conda:
        "../envs/rseqc.yaml"
    run:
        samples = whatever_loads_samples(input.samples)
        read_length = samples.loc[wildcards.sample, 'read_length']
        shell("deletion_profile.py -i {input.bam} -l {read_length} -o {output}")

```

(On a side note, you can use `{output}` rather than repeating yourself for that part.)

---
Tags: snakemake
---


title: Snakemake and github actions

text:
# Question
Title: Snakemake and github actions
I am working on a snake-make workflow that follows the guidelines of formatting and linting so that it can be published on the official Snake-make repository. Whenever I push my code on GitHub, the GitHub actions are run and the pipeline is tested as indicated here.

The problem is related to the execution of some time-consuming rules, such as the genome indexing generation for bowtie. These rules are run every time and take plenty of time. I was wondering whether there was a way to bypass the execution of these rules and test the pipeline anyway. Do you have any suggestions about this?

A problem is that I have a rule in my snakemake workflow that downloads the human genome and that is needed to index and then for the alignment. I need a workaround to avoid this behaviour, since has you said it is not feasible to test on the whole genome.

Thanks

# Answer
> 2 votes
The answer for any testing purpose is generally to put together a small and minimal dataset so the test completes in seconds or minutes while still being representative. It makes (to me) limited sense to have a Actions/CI running using for example the entire human genome as input. The point of testing is to check that steps are executed properly, arguments are passed correctly, errors are detected. For this use either a small genome (such as from a bacterium or yeast), or make an artificial one consisting of very few basepairs.

What I often do (for example for a peak calling workflow) is to take a full dataset, then subset to a region where I see a) many peaks and b) no peaks, then extract this genomic range (a few thousand bp) to build an artificial genome, plus extract from the bam/fastq the reads contributing to this region. That way I have basically a read-world dataset but very limited in size. For CI testing that does very well and CI finished in seconds, allowing to test even different configurations of this workflow on a standard / unpaid Actions runner.

I personally prefer to not automate downloads but rather to have the user provide an actual file, or provide the file pre-downloaded in the workflow bundle. Connections might be bad or broken that day or the server/HPC is not connected to internet, and the entire workflows fails. Maybe write the relevant section in a way that it either accepts a URL for download OR a local file, and then write a testing profile that then reads the file to ensure it always is accessable to the Actions runner.

---
Tags: python, snakemake, github
---


# Question
Title: Question about umap using different numbers of pca components as initialization
I am new to the scRNA-seq field and I have been doing some experiments of visualization of UMAP using different numbers of PCA components for initialization. The process involves projecting scRNA-seq data (count matrix) onto various numbers of PCA components, followed by non-linear dimension reduction using UMAP. It seems that the UMAP results are exactly the same across different initializations, regardless of the number of PCA components used. Are there any theoretical reasons for this? I am using Dynamo and their function `dyn.tl.reduceDimension(adata, n_pca_components=30)` for the experiments.

The results can be reproduced using the following code. Here I used the build-in zebrafish dataset with `n_obs × n_vars = 4181 × 16940`. I am just trying to see how the PCA initialization affects the Umap results.

```
import warnings
warnings.filterwarnings('ignore')

import dynamo as dyn

dyn.configuration.set_figure_params('dynamo', background='white')

adata = dyn.sample_data.zebrafish()

preprocessor = dyn.pp.Preprocessor(cell_cycle_score_enable=True)
preprocessor.preprocess_adata(adata, recipe='monocle')

dyn.pl.variance_explained(adata)

dyn.tl.reduceDimension(adata, n_pca_components=1, enforce=True)
dyn.pl.umap(adata, color='Cell_type')

dyn.tl.reduceDimension(adata, n_pca_components=30, enforce=True)
dyn.pl.umap(adata, color='Cell_type')

```

Here is the picture that shows the variance carried by PCA:

And the final umap results using n\_pca\_components= 1 and n\_pca\_components= 30:

# Answer
It seems you are performing PCA -\> UMAP. The results you are reporting actually look cool. It's a good answer and there's nothing to worry about. It also shows good analytics, 'cause you're questioning the robustness of your data.

The limiting factor is the variance carried by PCA - which needs to be calculated (it may be low, e.g. \<50%). Well I mean you've calculated it via the PCs - but its not declared, its better to know how much it is, i.e. how much variance each PC holds. The UMAP is doing nothing more than representing the established variance identified by PCA.

If you change the number of components which UMAP is assessing *then* this will result in changes in the clusters observed, usually.

**However** if you are changing the no. PCs -\> UMAP and there's little effect then the PCs are not introducing new information to cut-up the data set. Thats a fairly robust answer: it's attained "stationarity". Very desirable output, just make sure there isn't a bug in your code, but it demonstrates there's no hidden or partial information that might give a different result. So the clusters you are seeing are a good representation of the data (thats good).

Its not a complete analysis - as hinted - but I'd be happy with that result.

---

Overall, there's very heavy dispersion of the variance between PCs. I've never had to teach this stuff, beyond a few line explanations. An educator experienced with this theory would give a more succinct explanation. However, the discussion below arrives at a common understanding between the OP and I. If there's demand I'll summarise the discussion as a singular paragraph: the final comments should be the introduction, i.e. what PCA is doing at basic level.

**Comments** (discussion)

**@Zach** even when using n\_pca\_components=1 and n\_pca\_components=30, the results remain unchanged. The first PCA component should only capture a part of the overall information. Is this correct?

**Response** Number 0 (on your graph) is the first principle component (16%). PC = 1 is the second PC (5%). Thus PCA = 1 means 21% of the variance. This far exceeds the PC 3 at 2%. Thats why it works, so all information is PC1 and PC2, nothing beyond that discriminates against the data.

**@Zach** So, n\_pca\_components = 3 could capture almost all the information and adding more components doesn't really make any difference to the results. But if we only consider n\_pca\_components = 1, why there are still no differences compared to n\_pca\_components = 30?

**Response** its easy .. because if you add all the PCs from PC 4 + ... + PC 30 it will be \<5% of the total variance. PC1 + PC2 + PC3 = 24% of the variance, so there is nothing to see beyond PC3, when the maximum is PC 30.

**@Zach** Why there are no differences between n\_pca\_components = 1 (using only PC1) and n\_pca\_components = 30 (using PC1+PC2+...+PC30)?

**Response** PC1 + PC2 = 24% and PC1 ... + PC30 \<30%. There is no information in the variance for between group clusters for PC3 to PC30, there is likely to be within group differences but you are not measuring those.

**@Zach** there is little difference between PC1 + PC2 = 21% and PC1+...+PC30. I mean if we only consider PC1 = 16% and PC1+...+PC30 \> 25%, do they have any differences between PC1 and PC1+...+PC30?

**Response** PCA works on individual elements it, it doesn't work on groups. The analysis of groups is post-hoc (usually via colour coding). PC=1 in your UMAP calculation is PC1 + PC2 (PC1 is automatic, its position 0 in your graph). . You are only looking at between group variance and in your data that is resolved essentially for PC1 + PC2 = 21%. The within group positioning will change so, yes there will be difference against PC1+PC2 vs PC1 + .... + PC 30 but we can't see them. Again PCA is individual elements, it's not using group information. So if you did sub-group analysis you might see changes PC1+PC2 vs PC1 + ... PC 30. So to put it simply your group delineation is too big to spot in any differences (but that only matters if it is biologically relevant).

**@Zach** Actually what I mean is PC0 = 16% and PC0+...+PC30 \> 25%. But your explanation makes sense to me and I appreciate your response.
> 0 votes

# Answer
It is a bit odd that different numbers of principal components are giving identical results, especially given that even the same initialisation settings with different UMAP epoch numbers will give different results.

I'm not familiar with the use of python for single cell data analysis, but it looks like the `adata` information might be being modified when running `reduceDimension`. If this is the case, then it's possible that the first PCA reduction is reducing the dataset as far as it will go, and the second PCA reduction doesn't do anything. This could be checked by changing the order of operations (from a fresh run) and seeing if the output is changed:

```
dyn.tl.reduceDimension(adata, n_pca_components=30, enforce=True)
dyn.pl.umap(adata, color='Cell_type')

dyn.tl.reduceDimension(adata, n_pca_components=1, enforce=True)
dyn.pl.umap(adata, color='Cell_type')

```

It's also possible that the `umap` command is doing its own reduction and ignoring the result of `dyn.tl.reduceDimension`. Maybe the UMAP plotting function doesn't actually do a UMAP reduction as well, and the UMAP result (or initial condition) is calculated earlier (e.g. in `preprocessor.preprocess_adata`).

It's also possible, as @M\__ suggests, that you've reached a stationary result after however many UMAP iterations are set as default. But... I find that unlikely given how variable UMAP has been for me. I'd want to do a lot more exploration before settling on that as an explanation.
> 0 votes

---
Tags: rna-seq, pca, umap, dimensionality-reduction
---


title: How to subset samples from a VCF file?

text:
# Question
Title: How to subset samples from a VCF file?
I have VCF files (SNPs & indels) for WGS on 100 samples, but I want to only use a specific subset of 10 of the samples. Is there a relatively easy way to pull out only the 10 samples, while still keeping all of the data for the entire genome?

I have a script that allows me to pull out regions of the whole genome for all 100 samples, so if I could do something similar but only put regions for the 10 samples that I want that would be ideal.

# Answer
> 16 votes
Bcftools has sample/individual filtering as an option for most of the commands. You can subset individuals by using the `-s` or `-S` option:

> `-s`, `--samples [^]LIST`
>
> Comma-separated list of samples to include or exclude if prefixed with "^". Note that in general tags such as INFO/AC, INFO/AN, etc are not updated to correspond to the subset samples. `bcftools view` is the exception where some tags will be updated (unless the `-I`, `--no-update` option is used; see `bcftools view` documentation). To use updated tags for the subset in another command one can pipe from view into that command. For example:
>
> `-S`, `--samples-file FILE`
>
> File of sample names to include or exclude if prefixed with "^". One sample per line. See also the note above for the `-s`, `--samples` option. The command `bcftools call` accepts an optional second column indicating ploidy (0, 1 or 2) or sex (as defined by `--ploidy`, for example "F" or "M"), and can parse also PED files. If the second column is not present, the sex "F" is assumed. With `bcftools call -C trio`, PED file is expected. File formats examples:
>
> ```
> sample1    1
> sample2    2
> sample3    2
>
> ```
>
> or
>
> ```
> sample1    M
> sample2    F
> sample3    F
>
> ```
>
> or a .ped file (here is shown a minimum working example, the first column is ignored and the last indicates sex: 1=male, 2=female):
>
> ```
> ignored daughterA fatherA motherA 2
> ignored sonB fatherB motherB 1
>
> ```

Example usage:

```
bcftools view -s sample1,sample2 file.vcf > filtered.vcf
bcftools view -S sample_file.txt file.vcf > filtered.vcf

```

See the bcftools manpage for more information.

# Answer
> 4 votes
In addition to the answer from @gringer there is a `bcftools` plugin called `split` that can do this, but gives you the added ability to output single-sample VCFs by specifying a filename for each sample.

```
$ bcftools +split

About: Split VCF by sample, creating single-sample VCFs.

Usage: bcftools +split [Options]
Plugin options:
   -e, --exclude EXPR              exclude sites for which the expression is true (applied on the outputs)
   -i, --include EXPR              include only sites for which the expression is true (applied on the outputs)
   -k, --keep-tags LIST            list of tags to keep. By default all tags are preserved
   -o, --output DIR                write output to the directory DIR
   -O, --output-type b|u|z|v       b: compressed BCF, u: uncompressed BCF, z: compressed VCF, v: uncompressed VCF [v]
   -r, --regions REGION            restrict to comma-separated list of regions
   -R, --regions-file FILE         restrict to regions listed in a file
   -S, --samples-file FILE         list of samples to keep with second (optional) column for basename of the new file
   -t, --targets REGION            similar to -r but streams rather than index-jumps
   -T, --targets-file FILE         similar to -R but streams rather than index-jumps
Examples:
   # Split a VCF file
   bcftools +split input.bcf -Ob -o dir

   # Exclude sites with missing or hom-ref genotypes
   bcftools +split input.bcf -Ob -o dir -i'GT="alt"'

   # Keep all INFO tags but only GT and PL in FORMAT
   bcftools +split input.bcf -Ob -o dir -k INFO,FMT/GT,PL

   # Keep all FORMAT tags but drop all INFO tags
   bcftools +split input.bcf -Ob -o dir -k FMT

```

So if you had the following samples file `samples.tsv`

```
sample1    sample1
sample2    sample2

```

You can run it and get the following

```
$ bcftools +split -S samples.tsv -o outdir in.vcf
$ ls
in.vcf sample1.vcf sample2.vcf

```

Without the second column, you would just get a single VCF with the two samples in it (as you would with `view`)

# Answer
> 2 votes
You can use the GATK's SelectVariants tool with the `-sn` flag.

E.g.

`gatk SelectVariants -V input.vcf -R reference.fasta -sn Sample_01 -out sample.vcf`

You may use the `-sn` flag several times so as to select several samples, or use it to point to a file containing a sample name on every line.

# Answer
> 0 votes
From the old VCF manual,

https://vcftools.sourceforge.net/man\_latest.html

and

https://manpages.debian.org/testing/vcftools/vcf-subset.1.en.html

`cat in.vcf | vcf-subset -r -t indels -e -c SAMPLE1 > out.vcf`

comma-separated list of variant types ( indels ) could be obtained.

---
Tags: vcf, wgs
---


# Question
Title: How should I use dorado basecaller to calculate translocation time? Can I change models' config files?
I am trying to use **dorado** basecaller to process ONT data. I need to make the move table and I am doing it with **--emit-move**.

```
dorado basecaller dna_r9.4.1_e8_sup@v3.3/ /pod5_files/ --recursive --emit-moves --modified-bases 5mCG_5hmCG -v > /bams/mod_stride1.bam

```

The main purpose of the move table for me is to calculate the number of ONT device measurements before basecaller emits one base (calculate **translocation time per base**). But the stride value in config.toml file is fixed to 5 (stride=5). I need to calculate the instantaneous translocation time per base (in seconds), because I am looking to find a measure similar to IPD values in PacBio, when we can get the time for sequencing each base.

I was wondering to **change the stride value in config.toml file** in the model folder, from stride=5 to stride=1, to directly infer the number of measurements without using Tombo for correcting and re-squiggling.

Unfortunately, even though this command works perfectly fine when the stride=5, when I set stride=1 it always gives me the following error: **"CUDA out of memory"**, even when I triple the size of RAM for the job.

My questions are:

1. Does setting stride=1 make sense for this purpose? Or do I have to do something else? (I know we have some K-mer models too)
2. Is this even possible to change the stride value on config.toml file?
3. How should I figure out how much RAM I need for the job?

# Answer
In pod5 files the number of samples prior to the adapter/sequence should be encoded in the metadata, and can be looked at using ONT's `pod5 inspect` tool:

```
$ pod5 inspect reads pod5_file.pod5

read_id,channel,well,pore_type,read_number,start_sample,end_reason,median_before,calibration_offset,calibration_scale,sample_count,byte_count,signal_compression_ratio
00445e58-3c58-4050-bacf-3411bb716cc3,908,1,not_set,100776,374223800,signal_positive,205.3,-240.0,0.1,65582,58623,0.447
00520473-4d3d-486b-86b5-f031c59f6591,220,1,not_set,7936,16135986,signal_positive,192.0,-233.0,0.1,167769,146495,0.437
...

```

The average translocation speed can be worked out by comparing the `sample_count` to the number of called bases.

More details here.

Unfortunately, due to the way the basecaller models work, I don't think there's any easy way to create a 1-sample lookup table for precise instantaneous translocation calculations.

If it is just an out of memory error, then it should work (but very slowly) if you set the basecaller batch size to 1 (i.e. `-b 1`). If even that doesn't work, something else is going on and it's unlikely to be a memory problem.
> 1 votes

---
Tags: nanopore
---


# Question
Title: What options are there for visualising DNA sequences?
I typically see several ways of visualising DNA sequences and I would like to know when (best) to use them?

For example, I sometimes see the letter codes GATC, sometimes I see colours, sometimes I see a representation of the double-helix. I think I saw a circular representation.

Additionally, what other options exist for visualising DNA and when should they be used? E.g. for papers, analysis in software, etc.

What do people tend to use when presenting DNA sequences?

# Answer
> 2 votes
DNA sequences have a massive range in length, so the approach of visualisation is dependent on the purpose of that visualisation.

I have been experimenting with visualising the repetitive structure of DNA in a way that represents repeats at both short range (tens of bases) and long range (millions of bases, or longer). You can see my attempts at that here. An associated presentation about DNA visualisation (with emphasis on repetitive sequences) can be found here.

---
Tags: genome, visualization, dna
---


# Question
Title: How to access yeast-cyc in a programmatic way to extract the GO terms?
I am a newbie in Bioinformatics, I want to extract the GO terms from yeast-cyc, which is items in the Go terms tab of this page https://yeast.biocyc.org/gene?orgid=YEAST&id=G3O-29622#tab=GO , not one page but one database like *Saccharomyces cerevisiae S288c*.

Is there an API or python lib can help us finish this task? I had found biocyc in pypi, but I did not found Go terms was involved in this API.

If the API is not convinient, can I download some dataset in which my required items is include?

# Answer
> 1 votes
You can get the GO terms for yeast from UniProt, e.g. https://uniprot.org/uniprotkb?query=(taxonomy\_id:559292).

Using the UniProt website you can generate an API endpoint URL for your query, e.g. here are the GO terms for all yeast genes: https://rest.uniprot.org/uniprotkb/stream?fields=accession%2Cid%2Cprotein\_name%2Cgene\_names%2Cgo&format=tsv&query=%28%28taxonomy\_id%3A559292%29%29

You can download this data using `wget` or `cURL`, e.g.

```
wget -O yeast_go_terms.tsv "https://rest.uniprot.org/uniprotkb/stream?fields=accession%2Cid%2Cprotein_name%2Cgene_names%2Cgo&format=tsv&query=%28%28taxonomy_id%3A559292%29%29"

```

---
Tags: python, genomics, public-dataset
---


# Question
Title: How to calculate logCPM across all samples?
Using edgeR for differential analysis between Tumor and Normal gave me differential expressed genes with logFC, logCPM, PValue and FDR.

From the details of glmTreat function I see that logCPM is average log2-counts per million, the average taken over all libraries.

The output table I got from glmTreat:

```
Geneid     logFC  unshrunk.logFC    logCPM  PValue       FDR
Gene1   4.985298235 5.573931983 1.236660486 1.55E-54    1.87E-50
Gene2   4.613438634 5.52126484  0.920343233 2.59E-53    1.56E-49
Gene3   5.250601296 5.356432653 2.751294034 2.48E-50    9.95E-47
Gene4   4.943049159 5.379741165 1.361393757 2.06E-45    6.19E-42
Gene5   6.100121754 6.117580436 5.774385315 2.95E-43    7.11E-40
Gene6   4.722697891 5.320461275 1.120685402 3.05E-42    6.11E-39
Gene7   5.246129853 5.497012001 1.902992053 2.78E-40    4.78E-37
Gene8   3.878773277 4.956636276 0.776208741 8.66E-39    1.30E-35
Gene9   4.441752496 4.94930499  1.132652682 4.37E-38    5.83E-35

```

To calculate logCPM given in the above table manually I did like this:

```
logCPM <- cpm(y, prior.count=2,  log=TRUE)

```

I got the logCPM values for all the genes. Then I calculated Average across all the samples for each gene. I got like below:

```
Geneid  Average
Gene1   0.686560246
Gene2   0.617115826
Gene3   1.075975225
Gene4   0.692050878
Gene5   1.277556065
Gene6   0.638358189
Gene7   0.689323163
Gene8   0.60700396
Gene9   0.662115092

```

Why there is so difference? Where I'm doing wrong?

# Answer
If you want the values by group then subset `y` to contain the samples of interest and then feed that to `aveLogCPM()`.
> 3 votes

# Answer
This post on the Bioconductor mailing list explains the difference: https://support.bioconductor.org/p/100463/
> 0 votes

---
Tags: r, rna-seq, bioconductor, rna, edger
---


# Question
Title: Ordering average of gene expression in heat map
I want to plot a heat map in header I have the average of gene expression validated in lab for each sample

```
dput(heat[1:5,1:5])

structure(c(3.05574677859363, -0.451189292623471, -0.398066832453666,
-0.456927931251448, -0.393418432279742, -0.270234667496526, -0.451189292623471,
-0.398066832453666, -0.456927931251448, -0.393418432279742, -0.270234667496526,
-0.451189292623471, -0.398066832453666, -0.456927931251448, -0.393418432279742,
-0.270234667496526, -0.451189292623471, -0.398066832453666, 2.0338091136055,
-0.393418432279742, -0.270234667496526, 1.76078923919148, 1.96896026174603,
1.97960246116808, -0.393418432279742), .Dim = c(5L, 5L), .Dimnames = list(
    c("hsa-miR-548b-5p", "hsa-miR-181c-5p", "hsa-miR-181d-5p",
    "hsa-miR-125a-3p", "hsa-miR-4286"), c("30566-001", "30566-002",
    "30566-003", "30566-004", "30566-005")))

dput(m[1:5,c(2,7)])

structure(list(HealthStatus = c("Healthy", "Healthy", "Healthy",
"Healthy", "Healthy"), Average = c(108.915677, 103.148873, 40.3359607,
88.4329528, 87.7756577)), row.names = c("30566-001", "30566-002",
"30566-003", "30566-004", "30566-005"), class = "data.frame")

pheatmap(heat, cluster_rows=T, show_rownames=T,
         cluster_cols=FALSE, annotation_col=m[,c(2,7)],
         annotation_colors=annot_colors)

```

How I can order `average` in a descending way?

# Answer
> 4 votes
If you order your matrix/data.frame you can sort your columns by decreasing "Average", e.g.

```
library(pheatmap)

heat <- structure(c(3.05574677859363, -0.451189292623471, -0.398066832453666,
                    -0.456927931251448, -0.393418432279742, -0.270234667496526, -0.451189292623471,
                    -0.398066832453666, -0.456927931251448, -0.393418432279742, -0.270234667496526,
                    -0.451189292623471, -0.398066832453666, -0.456927931251448, -0.393418432279742,
                    -0.270234667496526, -0.451189292623471, -0.398066832453666, 2.0338091136055,
                    -0.393418432279742, -0.270234667496526, 1.76078923919148, 1.96896026174603,
                    1.97960246116808, -0.393418432279742), .Dim = c(5L, 5L), .Dimnames = list(
                      c("hsa-miR-548b-5p", "hsa-miR-181c-5p", "hsa-miR-181d-5p",
                        "hsa-miR-125a-3p", "hsa-miR-4286"), c("30566-001", "30566-002",
                                                              "30566-003", "30566-004", "30566-005")))

m <- structure(list(HealthStatus = c("Healthy", "Healthy", "Healthy",
                                     "Healthy", "Healthy"), Average = c(108.915677, 103.148873, 40.3359607,
                                                                        88.4329528, 87.7756577)), row.names = c("30566-001", "30566-002",
                                                                                "30566-003", "30566-004", "30566-005"), class = "data.frame")

pheatmap(heat[,order(m[,"Average"], decreasing = TRUE)],
         cluster_rows=TRUE, show_rownames=T,
         cluster_cols=FALSE, annotation_col=m)

```

Created on 2023-09-06 with reprex v2.0.2

---
Tags: heatmap
---


# Question
Title: Separate multiple sequence into different file, one sequence per file
I have a file contain multiple sequence, and I want to separate them by "gene:" into different file. example:

`example.fa`

```
>KQK21959
ATGGTGCTCGACCCCGTCGTCGGGGGCTCGGCCTCCCGCAAGGTCGAGCCCCTAACCCTA

>KQK14546
ATGCCTCCGCCGGACGCGGCCACCTACGCTCCCACGCTGCCAATGGCGACGTACCTCCCA

```

How to get results like this? One sequence per file, and after `>` is gene id, not cds id

`KQK21959.fa`

```
>KQK21959
ATGGTGCTCGACCCCGTCGTCGGGGGCTCGGCCTCCCGCAAGGTCGAGCCCCTAACCCTA

```

`KQK14546.fa`

```
>KQK14546
ATGCCTCCGCCGGACGCGGCCACCTACGCTCCCACGCTGCCAATGGCGACGTACCTCCCA

```

# Answer
A simple Biopython solution:

```
from Bio import SeqIO

for record in SeqIO.parse("example.fa", "fasta"):
    SeqIO.write(record, record.id + ".fa", "fasta")

```

Note, if you have multiple sequences with the same id, only the last sequence will be kept in the final output files
> 3 votes

# Answer
There's an app for that! The `exonerate` package (see here) includes a nice `fastaexplode` command that does exactly what you're looking for and will be much, much faster than anything we give you here since it is written in C and not a scripting language. Just install `exonerate` (you can do that with `sudo apt install exonerate` on Debian-based systems) and then:

```
fastaexplode file.fa

```

This will create these two files (`file.fa` is the original input file from your question):

```
$ ls
file.fa  KQK14546.fa  KQK21959.fa

```
> 3 votes

# Answer
That really depends on the tools you are using. If anything will do you can use the following python snippet.

If you know that each sequence is on a single line and your file isnt gigantic use this

# EDIT

Updated so the file isn't read to memory at once to address @terdon point. Although, the above solutions are probably much better.

```
with open( "seqs.fasta" ) as f:
    out = ""

    for line in f:
        if not line.strip(): continue
        if line.startswith( ">" ):
            if not type( out ) == str: out.close()
            gene = line[1:].strip().split()[0]
            out = open( f"{gene}.fa", "w" )
            out.write( line )
            continue
        out.write( line.strip() )

```
> 2 votes

# Answer
This can also be done using awk. This is harder to read, but probably faster than python:

```
awk 'sub(/^>/,""){
        FILE=$0".fa";
    print ">"$0 >FILE;
        next
    }{print > FILE}' file.fa.

```

This will define a new output file `FILE` at each header line, whose filename is the header minus the `>` character, suffixed with `.fa`. And all lines in the file are appended to the current output file.

EDIT: Changed to a faster version found on: https://crashcourse.housegordon.org/split-fasta-files.html#using-awk-instead-of-csplitsed
> 2 votes

# Answer
This worked for me today:

```
# Load required packages
BiocManager::install("Biostrings")

library(Biostrings)

# Read the multi-genome FASTA file
multi_genome_fasta <- readDNAStringSet("*.fasta")

# Extract individual genomes and save as separate FASTA files
for (i in 1:length(multi_genome_fasta)) {
    header <- names(multi_genome_fasta)[i]
    sequence <- as.character(multi_genome_fasta[i])

    # Combine header and sequence
    combined_data <- paste(">", header, sequence)

    # Save the combined data as a separate FASTA file
    file_name <- paste0(header, ".fasta")
    writeLines(combined_data, file_name)

    cat("Saved", header, "as", file_name, "\n")
}

cat("Separation completed.\n")

```
> 1 votes

# Answer
There are many ways to do it. Here is one way to do in Python.

```
import os

with open("input.fasta") as f:
    for line in f:
        if line.startswith('>'):
            record_name = line.lstrip(">")
            outfilename = f"{(record_name.rstrip())}.fasta"
            if not os.path.exists(outfilename):
                out = open(outfilename, 'w+')
            else:
                print('File already exists!')
        out.write(line)

```
> 1 votes

---
Tags: fasta, ngs, gene
---


# Question
Title: Why does the subset of rows that match a condition added by the rows that don't match a condition not add up to the total number of rows?
I need to match the external gene names of a table I have called (which has 1,000,000 entries) with their respective ensemble gene IDs.

Here's a small subset of `df`:

| V1 | external\_gene\_name | V3 |
| --- | --- | --- |
| M0085 | A1BG | 1 |
| M0085 | SNX8 | 1 |

and here's a small subset of `ensembl_ids`:

| ensembl\_gene\_id | external\_gene\_name |
| --- | --- |
| ENSG00000288263 | GRK1 |
| ENSG00000291428 | BRSK2 |

I also included a line to see the number of entries in `external_gene_name` which do not have corresponding ensemble names. Here is my code:

```
merged_df <- merge(df, ensembl_ids, by = "external_gene_name")
elements_not_in <- df$external_gene_name[!(df$external_gene_name %in% ensembl_ids$external_gene_name)]

```

But what's confusing me is that when I add the number of rows in `merged_df` with that of `elements_not_in`, the total exceeds 1,000,000. Shouldn't the total be 1,000,000? I'm wondering if there's a flaw in my code.

# Answer
> 4 votes
From a brief glance, the gap is not in your code but in your expectations. You're operating based on an assumption of 1:1 mapping between ENSG IDs and HGNC symbols (your `external_gene_name` attribute). In reality, multiple ENSG IDs will map to the same HGNC symbol so the merge will produce multiple rows per `external_gene_name` in the merge product. This is why you're seeing a total greater than the sum of number of unique entries when you compare the merge product + the anti-join result (i.e., the `elements_not_in` dataset) to the initial individual datasets.

---
Tags: r, ensembl
---


# Question
Title: How to analyze GC content of fasta file?
I have the assembly file of a certain lizard species, I wanted to recreate Fig. 2 of the following paper, which compares the GC content distribution over 5kb windows, how would I go about doing this?

https://bmcecolevol.biomedcentral.com/articles/10.1186/s12862-022-02086-7#Sec10

# Answer
> 4 votes
I would suggest using `bedtools makewindows` and `bedtools nuc` to compute the GC content for genomic windows as shown in Fig 2 of the linked paper.

See related question here.

This will slice the genome into 5Kbp (or whatever size specified to `makewindows`) pieces and compute the GC content of each piece, which can be plotted as a distribution as in the figure.

Something like this would work for the 5000 bp example:

```
# assemble a bedtools compatible genome file:
# see https://www.biostars.org/p/70795/ for hints how to do this.
bedtools makewindows -w 5000 -g my_genome.genome > windows.bed # you just made this .genome file manually
bedtools nuc -bed windows.bed -fi my_genome.fasta > windows_gc.bed

```

The bed output file will contain the GC content of each window in the `7_pct_gc` column (according to my understanding of the output file).

---
Tags: fasta, genome
---


# Question
Title: How to do average macrophages signature boxplots, with Bulk RNA data from 2 control and 2 treatment samples and LM22 signature matrix from cibersort?
I have bulk data from 4 samples of triple negative breast cancer, which 2 of them are control and the other 2 are treated. The goal here is to estimate the average macrophages signatures, for M0,M1 and M2 in each of the 4 samples.

My bulk data:

```
Gene    MCtrl1_TPM  MCtrl2_TPM  MTreat1_TPM MTreat2_TPM
KHDC1L  0.0 75.15   0.0 0.0
CD24    1.53    41.56   0.11    0.06
MAMLD1  5.2 11.94   0.0 0.0
SMIM31  4.1 10.24   0.0 0.06
C16orf95    7.9 13.24   0.0 0.08
OCLN    2.69    7.06    0.06    0.03
RBM14-RBM4  16.97   1.95    0.0 0.0
RPS10-NUDT3 0.0 0.0 31.91   12.06
NDUFC2-KCTD14   15.06   49.98   0.41    0.0
ATP9A   2.7 10.05   0.01    0.0
UST 0.38    10.63   0.0 0.01
LOC100996709    12.25   30.58   0.0 0.0
ADAM8   10.13   2.14    0.03    0.0
ADGRG2  0.05    21.49   0.0 0.0
MBNL2   4.61    5.04    0.0 0.0
PIGBOS1 56.73   51.33   0.6 0.09

```

To do that, I thought was to use the CIBERSORTx software more specifically, the option - Impute Gene Expression, High-Resolution Mode (tutorial 5) with my bulk data and signature matrix file, the LM22 signature matrix file of 22 immune cells provided by CIBERSORTx, but with the selected columns correspondent to the 3 types of macrophages.

Signature matrix file:

```
Gene    M0  M1  M2
ABCB4   7.951360487 27.55771004 121.4322774
ABCB9   56.25188749 29.38087813 40.39026017
ACAP1   48.97221651 30.97425375 23.05515036
ACHE    13.6929629  1373.06962  48.58576071
ACP5    13644.1723  1062.475909 7633.960182
ADAM28  54.96648703 274.639997  360.6820442
ADAMDEC1    4464.262133 8524.022801 5004.969419
ADAMTS3 45.56049098 32.34016681 33.35142773
ADRB2   143.0178779 787.6239741 363.1850058
AIF1    1841.269508 1611.228571 11610.29035
AIM2    181.9571    2965.888162 148.2845223
ALOX15  6.621226619 5.412216318 430.3353048
ALOX5   788.5397751 83.30089783 248.3724724
AMPD1   76.54456664 62.39764428 59.95269372
ANGPT4  39.33257003 21.218127   11.8097012
ANKRD55 123.600429  96.20392388 48.55417359
APOBEC3A    85.2400852  4501.140479 148.7993403
APOBEC3G    278.7681406 1036.993881 258.3004281
APOL3   238.5980956 11009.29392 437.5272083
APOL6   13.10943681 967.6451456 7.784612131
AQP9    4583.500675 4309.341662 399.4746192
ARHGAP22    175.4574354 351.5922058 272.9890572
ARRB1   333.3962724 946.0138184 395.5562302
ASGR1   108.7218916 25.23626149 155.4913385
ASGR2   36.11978652 26.3738691  760.7257707
ATHL1   52.25820981 167.5032595 132.631945

```

But I encountered this error, trying to do that:

```
## Formatted for easier reading
Error in runCIBERSORTxGEP(weights, mixture, classes, QN, maxsamples, label, :
Abort! The number of mixture samples is not greater than the number of cell types.
Calls: CIBERSORTxHiRes -> CIBERSORTxGEP -> runCIBERSORTxGEP Execution halted

```

How this error is possible, when I have 3 types of macrophages and 4 samples? There's another way to estimate the average of the 3 types of macrophages signatures? Specifically in a graphic like this?If yes, how?

# Answer
> 2 votes
You can use `IOBR` R package for this purpose. `IOBR` is an R package to perform comprehensive analysis of tumor microenvironment and signatures for immuno-oncology.

```
if (!requireNamespace("remotes", quietly = TRUE))
install.packages("remotes")

if (!requireNamespace("IOBR", quietly = TRUE))
remotes::install_github("IOBR/IOBR")

library(IOBR)
library(tidyverse)

file = read.delim(file.choose(),row.names=1)        #bulk data from 4 samples

# estimating the average macrophages signatures using cibersort method
cibersort<-deconvo_tme(eset = file, method = "cibersort", arrays = TRUE, perm = 100 )

head(cibersort)

#barplot
res<-cell_bar_plot(input = cibersort[,c(1,15:17)], title = "CIBERSORT Cell Fraction")

```

---
Tags: r, cibersortx
---


# Question
Title: missed cleavages combinations generation
I am looking for an R based solution to the following problem:

Trypsin does not cleave perfectly in proteomics experiments. So starting from a protein sequence FASTA I want to generate all possible peptides which COULD be generated if trypsin was allowed n number of "missed cleavages"

In informatics terms I would describe this as skipping a split during a `strsplit()` function

Example to clarify:

We start out with the protein sequence: `GRGKA` which is usually split by `strsplit()` (digested by trypsin) after every K and R

```
"GRGKA" %>% strsplit(., "(?<=K|R)", perl = TRUE)

[[1]]
[1] "GR" "GK" "A"

```

now imagine that trypsin did not cleave perfectly, for example missing always one K or R

```
"GRGK" "GKA"

```

would be generated in this case. Or if trypsin digested this protein really badly and would end up missing two cleavages we would end up with the sequence:

```
"GRGKA"

```

I would love to have a generalizable solution that allows to generate all combinations of sequences that arise with n number of missed cleavages.

If you need any more information please feel free to ask, sorry if I missed something!

jay.sf's requested example:

```
"FLGERTCTRNERQP"

```

perfect digest:

```
"FLGER" "TCTR" "NER" "QP"

```

one missed cleavage all combinations:

```
"FLGERTCTR" "TCTRNER" "NERQP"

```

two missed cleavages all combinations:

```
"FLGERTCTRNER" "QP" "FLGER" "TCTRNERQP"

```

three missed cleavages (there are a total of three K and Rs so trypsin would miss all)

```
"FLGERTCTRNERQP"

```

For longer sequences I think this becomes a frameshift like problem: If we count the K and Rs in the sequence:

FLGE**R**TCT**R**NE**R**QPFLGE**R**TCT**R**NE**R**QP
FLGE**1**TCT**2**NE**3**QPFLGE**4**TCT**5**NE**6**QP
FLGE**R**\_TCT**R**\_NE**R**\_QPFLGE**R**\_TCT**R**\_NE**R**\_QP (perfect digest)
FLGE**R**TCT**R**\_NE**R**QPFLGE**R**\_TCT**R**NE**R**\_QP (one missed cleavage and always the uneven R and Ks)
FLGE**R**\_TCT**R**NE**R**\_QPFLGE**R**TCT**R**\_NE**R**QP (one missed cleavage and always the even R and Ks)

So I think for 1 missed cleavage we have two frames over the whole sequence and for two missed cleavages we would have three frames and so on....

# Answer
Essentially we want to combine *1, 2, ..., n - 1* splits, and if the `var`iance of split lengths is greater than one we have a second solution. Here is a function that, after the `strsplit`, first creates indices `cb` to later `paste` the respective splits together. I'm not quite sure what the maximum number of splits is, but this might be a start.

```
mccg <- function(x) {
  ## missed cleavages combinations generation ++++++++++++++++++++++++++++++++++
  gs <- \(n, k, length=n*k, right=FALSE) {
    ## helper, similar to `gl()` but w/o factorization
    if (right) {
      cbr <- gs(n, k, n, right=FALSE)
      rev(max(cbr) - cbr + 1L)
    } else {
      rep_len(rep.int(seq_len(n), rep.int(k, n)), length)
    }
  }
  ## split string after [K, R]
  spl_x <- el(strsplit(x, '(?<=K|R)', perl=TRUE))
  len <- length(spl_x)
  res <- if (len <= 1L) {
    ## if [K, R] at the end or missing
    x
  } else {
    ## nested list of combinations
    cb <- lapply(1:(len - 1L) |> {\(.) setNames(., . - 1L)}(), \(i, .len=len) {
      cb <- list(gs(.len, i, .len))
      if (i != 1L && var(tabulate(el(cb)))) {
        ## second solution, starting from right
        c(cb, list(gs(.len, i, .len, right=TRUE)))
      } else {
        cb
      }
    }) |> c(setNames(list(rep_len(1L, len)), len - 1L))
    ## paste according to sequences
    rapply(cb, \(x) vapply(unname(split(spl_x, x)), paste, collapse='',
                           character(1)),
           how='list')
  }
  return(res)
}

```

## Usage

*Two splits:*

```
mccg(x='GRGKA')
# $`0`
# [1] "GR" "GK" "A"
#
# $`1`
# $`1`[[1]]
# [1] "GRGK" "A"
#
# $`1`[[2]]
# [1] "GR"  "GKA"
#
#
# $`2`
# [1] "GRGKA"

```

*Three splits:*

```
mccg(x='FLGERTCTRNERQP')
# $`0`
# [1] "FLGER" "TCTR"  "NER"   "QP"
#
# $`1`
# [1] "FLGERTCTR" "NERQP"
#
# $`2`
# $`2`[[1]]
# [1] "FLGERTCTRNER" "QP"
#
# $`2`[[2]]
# [1] "FLGER"     "TCTRNERQP"
#
#
# $`3`
# [1] "FLGERTCTRNERQP"

```

*More splits:*

```
mccg(x='FLGERTCTRNERQPFLGERTCTRNERQP')
# $`0`
# [1] "FLGER"   "TCTR"    "NER"     "QPFLGER" "TCTR"    "NER"     "QP"
#
# $`1`
# $`1`[[1]]
# [1] "FLGERTCTR"  "NERQPFLGER" "TCTRNER"    "QP"
#
# $`1`[[2]]
# [1] "FLGER"       "TCTRNER"     "QPFLGERTCTR" "NERQP"
#
#
# $`2`
# $`2`[[1]]
# [1] "FLGERTCTRNER"   "QPFLGERTCTRNER" "QP"
#
# $`2`[[2]]
# [1] "FLGER"          "TCTRNERQPFLGER" "TCTRNERQP"
#
#
# $`3`
# $`3`[[1]]
# [1] "FLGERTCTRNERQPFLGER" "TCTRNERQP"
#
# $`3`[[2]]
# [1] "FLGERTCTRNER"     "QPFLGERTCTRNERQP"
#
#
# $`4`
# $`4`[[1]]
# [1] "FLGERTCTRNERQPFLGERTCTR" "NERQP"
#
# $`4`[[2]]
# [1] "FLGERTCTR"           "NERQPFLGERTCTRNERQP"
#
#
# $`5`
# $`5`[[1]]
# [1] "FLGERTCTRNERQPFLGERTCTRNER" "QP"
#
# $`5`[[2]]
# [1] "FLGER"                   "TCTRNERQPFLGERTCTRNERQP"
#
#
# $`6`
# [1] "FLGERTCTRNERQPFLGERTCTRNERQP"

```

## Simulation

Here are some toy data,

```
set.seed(42)
n <- 1e3
min_chr <- 5; max_chr <- 50
strings <- replicate(n, {
  repeat {
    smp <- sample(setdiff(LETTERS, c("B", "Z", "J", "O", "U", "X")),
                  sample(min_chr:max_chr, 1L), replace=TRUE)
    if (any(grepl('K|R', smp))) break
  }
  paste(smp, collapse='')
}) |> setNames(nm=_)

```

over which the execution of the function can be tested in a loop.

```
r <- lapply(strings, \(x) tryCatch(mccg(x), error=\(e) NA))

```
> 2 votes

# Answer
**1)** gsubfn is like gsub except the replacement can be a proto object containing a `fun` proto method which is a function whose first argument receives the proto object and which inherits a `count` (and other) built-ins. The match is input to the function as the second argument and is replaced by the output of the function. Here we replace each R or K with itself and if it is more than the k-th occurrence then with a semicolon as well. Then split by semicolon.

```
library(gsubfn)

# inputs
x <- c("GRGKA", "GRGKAK")  # vector of inputs
k <- 1 # no of splits to skip

p <- proto(fun = function(self, x) paste0(x, if (count > k) ";"))
strsplit(gsubfn("[RK]", p, x), ";")

```

giving

```
[[1]]
[1] "GRGK" "A"

[[2]]
[1] "GRGK" "AK"

```

**2)** Insert a semicolon after each R or K and then remove the first k. `strsplit` what is left. This uses base R only.

```
xx <- gsub("([RK])", "\\1;", x)
for(i in seq_len(k)) xx <- sub(";", "", xx)
strsplit(xx, ";")
# same output as above

```
> 0 votes

---
Tags: r, proteins, sequence-analysis
---


# Question
Title: Error output in nextflow pipeline using fasterq-dump
I have one file with multiple SRA accession called `Dataframe_with_accession.txt`, for this example I just put one SRA in that file called :`SRR8933535`

And the idea is to create a nextflow pipeline to download the sra files in the `Dataframe_with_accession.txt`, zip the files using `pigz` and remove the previous `fastq file`

To do that I used this code :

```
params.inputFile = "/crex/proj/Dataframe_with_accession.txt"
params.outputDir = "/crex/proj/Output"

process DOWNLOAD_FASTQ {

    publishDir params.outputDir, mode: 'symlink'

    input:
    each (accession)

    output:
    tuple path("params.outputDir/${accession}_1.fastq.gz")
    path("params.outputDir/${accession}_2.fastq.gz")

    script:
    """
    fasterq-dump --threads 12 --outdir ${params.outputDir}  ${accession}
    pigz -p12 ${params.outputDir}/${accession}_1.fastq
    pigz -p12 ${params.outputDir}/${accession}_2.fastq
    """
}

workflow {
    // Read accessions from the input file
    accList = file(params.inputFile).readLines()

    DOWNLOAD_FASTQ(accList)
}

```

Then I run this nextflow file.

But I got the following error message :

```
[01/298f4b] process > DOWNLOAD_FASTQ (1) [100%] 1 of 1, failed: 1 ✘
Error executing process > 'DOWNLOAD_FASTQ (1)'

Caused by:
  Missing output file(s) `SRR8933535_1.fastq` expected by process `DOWNLOAD_FASTQ (1)`

Command executed:

  fasterq-dump --threads 12 --outdir /crex/proj/Output  SRR8933535

Command exit status:
  0

Command output:
  (empty)

Command error:
  spots read      : 8,389,148
  reads read      : 16,778,296
  reads written   : 16,778,296

Work dir:
  /crex/proj/Output/work/01/298f4be6dde2ee1856df31db169bec

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

```

And I do not undestrand what is going wrong here ? It says

```
Command output:
  (empty)

```

But the file is created

# Answer
> 2 votes
Each nextflow process looks in its working directory for the files listed in `output:` . In your `script:` you already direct the output to its eventual destination. After running the script nextflow tries to find the output files and stops, because it can't find them. Solution: keep the output files in your workdir, nextflow will detect them and use the `publishDir` directive to copy them to their eventual destination.

---
Tags: nextflow, fastq-dump
---


title: Blastp MSA to the same length

text:
# Question
Title: Blastp MSA to the same length
To generate aligned protein sequence truncate to the same length, I use blastp in biopython for MSA. I have multiple queries and one subject sequence. However, the alignment never returns the same length.

How can I align all sequences altogether, thus resulting in a set of sequences of the same length?

# Answer
> -1 votes
@Matteo Ferla wrote in the comments,

> It sounds like you'd want a MSA of the whole set and not a set of pairwise ones. If global MSA is what you actually want, try ebi.ac.uk/Tools/msa/muscle. Then to trim the ends before phylogeny there are dozens of tools.

> The reason is because Blast aligns locally, i.e. a window, while a global alignment is the whole entirety of the sequence

# Answer
> 2 votes
BLAST isn't appropriate for multiple sequence alignment. BLAST generates pairwise alignments. Furthermore, the aligning regions shown in a BLAST report only include subsets of the hit sequences.

To create multiple sequence alignment of hits you want to:

* run blast
* extract the hit sequences
* perform multiple-sequence-alignment (e.g., using MAFFT)

Extracting hit sequences can be done using the command-line. Check the docs for `blastdbcmd`.

If you prefer a point-and-click interface, you can also do this using my lab's SequenceServer BLAST software by simply selecting the hits that you want to retain, and in the **Download** section, clicking "FASTA of selected hits". Or if you want all of them, click "FASTA of all hits".

---
Tags: python, blastp
---


# Question
Title: Predict population based on PC coordinates
Utilizing a reference panel, I want to assign most likely population label to each individual in the study. Following are the files I have:

1. Reference panel population labels:

```
HG00121 EUR
HG00122 EUR
HG00123 EUR

```
2. Reference panel PC (Principal Component) coordinates:

```
HG00121 -0.0108487      0.0259052       -0.0099605      0.0168545
HG00122 -0.0115789      0.0264642       -0.01722        0.0157885
HG00123 -0.00996687     0.0266983       -0.0155755      0.0155963

```
3. Study individual PC coordinates:

```
Study_individual1 -0.0108487      0.0259052       -0.0099605      0.0168545

```

How do I go about assigning a population label to the study individual based of the reference panel? Most of the softwares I have come across such as evaluate\_check\_ancestry(plinkQC), FRAPOSA rely on inputs provided to PLINK which is not helpful in my case.

# Answer
> 1 votes
I apologise for being negative about this, but I am afraid you cannot use principle components for population inference other than assigning clusters within a bivariate or 3D graph. Eigen values are complex derivatives of data to achieve dimensionality reduction.

For population genetics you require a sample size =\> 10. I don't know the organism nor the question but goal is to derive allelic frequencies for each genome position. A common metric is F-statistics, for example Fst or Hardy-Weinberg equilibrium.

# Answer
> 1 votes
Here is a very simple way to get your started. Of course you can make it more complex:

Take the mean PC1 and PC2 coordinates for each 'reference panel population label' and then calculate the reference panel population which has the closest Cartesian distance to your target, i.e. $\sqrt{{(PC1\_{r} - PC1\_{t})^2 + (PC2\_{r} - PC2\_{t})^2}}$, where $PC1\_{r}$ is the mean PC1 for a ref population and $PC1\_{t}$ is PC1 for the target.

---
Tags: pca, population-genetics
---


# Question
Title: paired-end short reads: will one file suffice?
I have a quick question about paired-end short reads. I have multiple genomes that were sequenced with paired-end Illumina NextSeq 200 technology, resulting in two fastq files per sample: `001_R1.fastq.gz` and `001_R2.fastq.gz`, typical for short read data.

I would like to use a k-mer counting and clustering approach to put these genomes together on a phylogenetic tree. I was able to use both reads in one go to count k-mers but, as you can probably guess, this resulted in a very high frequency of k-mers identified and accounted for.

For this specific use case (not interested in changing my approach), would it make more sense to just use a single file for kmer counting (i.e. only the `001_R1.fastq.gz` file)?

**EDIT:** For more context, I'm trying to do an alignment-free approach to produce a phylogenetic distance comparison of these strains. I was using long reads before and it worked fine but Nanopore reads tend to be very messy. I have short-read data for all of the strains and would like to use that instead but that comes at the price of an overload of data.

I suppose my main concern/question is how do I pre-process these short reads without assembling or aligning to a reference to count kmers? Or is that not possible without at least alignment?

# Answer
> 3 votes
I suspect you won't be able to avoid assembling or aligning. If you are working with strains that don't have their own, well studied reference genome and so want to avoid aligning so as not to bias your data towards whichever strain has been chosen as the reference for your species (please correct me if this isn't what's going on), then you will need to assemble a denovo genome from these reads. At the very least you will need to build the longest consensus contigs you can. Unfortunately, this isn't something I have experience with so I can't help you do this beyond telling you that, as I'm sure you already know, there are tools out there that will try to assemble your reads into contigs.

I will instead try to explain why I think you can't avoid assembly. Short read fastq files, paired or not, contain short sequences of ~150-200 base pairs each. Now, depending on the depth of coverage, each region of your genome will be covered by multiple such short reads. Therefore, the data you have in your fastq files will be extremely repetitive and redundant since each of your genome's regions will be represented by multiple reads, and any attempt to find k-mers will be dogged by this redundancy. This is why you got so many more results than you expected and I don't see how using one file only would help at all.

If you had sequenced using an assay that includes unique molecular identifiers (UMIs), then you would be able to extract non-redundant, representative reads for all regions but I cannot think of a way to do this without UMIs off the top of my head. Not saying it isn't possible, only that you need someone more clever than I.

That said, the best I can come up with is that **if** your samples were sequenced in the same sequencing run and so can be assumed to have the same sequencing biases, you might be able to build your trees despite the redundancy in the data. My thinking here is that what you are really interested in isn't the k-mers themselves or their level of over-representation in each species, but only the difference in the k-mer content across your samples. If so, given that you sequenced all samples in the same run, they should all have the same level of over-representation and any differences might be enough for you to be able to cluster/build your trees.

This seems worth a try to me, but I haven't worked on something like this in many years so if you go down this route, make sure to double check and ensure I'm not leading you astray!

# Answer
> 1 votes
The problem with using single-end short reads is that they're often not long enough to uniquely match to a database (especially for highly-conserved genes, like the ribosomal RNAs). Whether or not this is a problem for your particular situation is unclear; you haven't provided enough context / information for me to work that out. For example, a single end of a 2x50bp run on a NextSeq 2000 might cause problems (especially for medium to large eukaryote genomes), but a 2x300bp run on a NextSeq 2000 is probably okay.

FWIW, Kraken2 has a `paired` mode, which combines the results from both read ends from the same template source:

> Usage of --paired also affects the --classified-out and --unclassified-out options; users should provide a # character in the filenames provided to those options, which will be replaced by kraken2 with "\_1" and "\_2" with mates spread across the two files appropriately. For example:

```
kraken2 --paired --classified-out cseqs#.fq seqs_1.fq seqs_2.fq

```

See the Wiki / manual here:

https://github.com/DerrickWood/kraken2/wiki/Manual#classification

---
Tags: sequence-analysis, illumina, k-mer
---


# Question
Title: Identifying the promoter on the reverse strand using visualization software (seqmonk)
I feel like this might be a really silly question and I've been going back and forth for awhile and I think confusing myself- so sorry if this is dumb.

I am using Seqmonk to visualize genes in an assembly (but I suspect the issue would be the same in any visualization software). I am trying to identify promoters of specific genes to determine if they are differentially methylated. Genes on the forward strand (highlighted in red), the promoter would inherently be upstream- in this case to the left (i.e. if my start site was 5000 and the end site 6000, I would assume I could find my promoter to the left of the start site between 3000 and 5000 bases.

My confusion is what happens with genes that are on the reverse strand (in blue here). Start site and end site could be the same (5000 and 6000- as end sites are always bigger than the start site). Does this mean that in the software, the promoter is expected to be on the right hand side (i.e. bases 6000-8000)? This has been my going thought but now I wonder if having bases that seem downstream (because they're bigger) makes sense, even though it really does because downstream on the forward strand would be upstream on the reverse strand. Or if I need to always do 3000-5000 bases as the promoter, regardless of the the strand of origin.

# Answer
> 2 votes
Have a look at this UCSC browser shot.

Looking at the *WDR92* transcript, note that it is transcribed right to left. This would indicate that the promoter is on the right. Indeed, we see a variety of promoter histone marks on the right of the gene, but not on the left.

In other words, yes, for genes on the reverse strand, the larger coordinate will be at the beginning of the gene, and the promoter will include a larger coordinate than the rightmost transcribed base.

---
Tags: genome, sequence-annotation, visualization
---


# Question
Title: How can I downsample a BAM file while keeping both reads in pairs?
I know how to downsample a BAM file to lower coverage. I know I can randomly select lines in SAM, but this procedure can't guarantee two reads in a pair are always sampled the same time. Is there a way to downsample BAM while keeping pairing information intact?

# Answer
> 17 votes
samtools has a subsampling option:

> -s FLOAT: Integer part is used to seed the random number generator \[0\]. Part after the decimal point sets the fraction of templates/pairs to subsample \[no subsampling\]

```
samtools view -bs 42.1 in.bam > subsampled.bam

```

will subsample 10 percent mapped reads with 42 as the seed for the random number generator.

# Answer
> 4 votes
With this function, you can subsample a given BAM file to a user-specified number of reads using SamBamba. The function automatically calculates the scaling factor. Ideas in part taken from here

```
function SubSample {

## Calculate the sampling factor based on the intended number of reads:
FACTOR=$(samtools idxstats $1 | cut -f3 | awk -v COUNT=$2 'BEGIN {total=0} {total += $1} END {print COUNT/total}')

if [[ $FACTOR > 1 ]]
  then
  echo '[ERROR]: Requested number of reads exceeds total read count in' $1 '-- exiting' && exit 1
fi

sambamba view -s $FACTOR -f bam -l 5 $1

}

## Usage example, selecting 100.000 reads:
SubSample in.bam 100000 > subsampled.bam

```

Note that `$2` is the number of reads. Say you choose 10000 reads but have paired-end data, you'll end up with 5000 fragments.

# Answer
> 1 votes
A totally different tactic...it helps if you aren't super picky about how many reads you have at the end.

Illumina read names contain the lane, tile, and xy coordinates of the read in them. If you filtered by these, you would be guaranteed to get either both members of a pair, or neither. I would recommend trying to avoid the edges of the flowcell.

# Answer
> 0 votes
It have used both sambamba and samtools. Either one works great to downsample to a desired read depth after calculating the downsample factor.

Sambamba can use multiple threads by using `-t INT`, where INT is the number of threads. It can also set the compression level for the BAM file using `-l`.

```
BAM="Sample.bam"
READS="10000000"
fraction=$(samtools idxstats $BAM | cut -f3 | awk -v ct=$myReads 'BEGIN {total=0} {total += $1} END {print ct/total}')

# With samtools
samtools view -b -s  ${fraction} ${BAM} > output.bam

# With sambamba
sambamba view -f bam -l 5 -s  ${fraction} ${BAM} > output.bam

```

---
Tags: sam
---


# Question
Title: DATES: sh: 1: dates_jackknife: not found
I've been trying to run DATES, a software that estimates admixture time, on my Linux Mint machine but I keep running into problems. This is the code to run the program:

```
./dates -p pardates.txt

```

The program then runs smoothly, but right before it ends it comes into an error:

```
Calling fit!...
dates_jackknife -p pardates.txt -l     0.450 -z Assyrian -r 77  -a >Assyrian:log
**sh: 1: dates_jackknife: not found**
##end of dates:       28.806 seconds cpu      627.270 Mbytes in use

```

I have tried the software on other Linux flavors but it keeps giving the same error, despite the fact I have exported the path to the bin directory every single time. I exported the path with this line:

```
export PATH=/ubuntu/home/Desktop/DATES/bin

```

I've also checked my version of perl script (dates jackknife's script) and it's up to date.

I have contacted the authors about the issue but they told me it's an issue with my environment, which I'm pretty sure is not.

# Answer
I suspect the path should not be 'ubuntu/home/Desktop/DATES/bin'. It should be /home/username/DATES/bin, where username would be your ubuntu ID/name
> 1 votes

---
Tags: genetics, admixture
---


# Question
Title: ETE4 manual or improved code description
ETE3 ***** (Environment for Tree Exploration) is a major package and in particular Python package to navigate phylogenetic trees. It is the premier coding library/package for tree manipulation.

ETE4.0.0 beta has been released via GitHub. The download contains an old ETE manual (2009) and a README file for ETE4, thats about it. Does anyone know the location of the new ETE4 manual, or list of functionality beyond what is in the README (posted below). Reading the code is the last resort ;-)

---

README ETE4

* Read and write support for trees in Newick format
* Multiple functions for traversing, searching and manipulating tree topologies and node annotations
* Integration with NCBI Taxonomic database
* Integration with GTDB database
* Programmatic visualization framework
* Visualization of huge trees based on adaptive zooming
* Comparing trees
* Phylogenetic functions:
    + orthology detection
    + phylogenetic distance?

---

**Reference** ***** ETE 3: Reconstruction, analysis and visualization of phylogenomic data. Jaime Huerta-Cepas, Francois Serra and Peer Bork. Mol Biol Evol 2016; doi: 10.1093/molbev/msw046

# Answer
> 1 votes
Just to feedback, the GitHub ticket was:

> ETE4.0.0 beta contains an old ETE manual (2009) and a README file for ETE4. Is there an updated ETE4 manual beyond etetoolkit.org/docs/latest/index.html or the README of ETE4?

lm-jkominek wrote,

> I know it's a bit of necro, but I wanted to second this request - while the current docs work fine enough for ETE3, they have diverged quite far from ETE4 and there's a bunch of errors that come through when using the latest version that have no explanation...

Thus I'm not alone ... (nothing to do with UAPs BTW)

jhcepas wrote,

> absolutely. We are finishing some final changes to the core of ETE4 and will start with the documentation. In the mean time, @jordibc has a list of most important changes regarding compatibility. Could you share somewhere as a wiki page @jordibc ?

jordibc replied,

> Yes, we are working on an updated version of the manual to be released with ete4. In the meantime, most of the old code can be run with some of the changes mentioned in the part of the wiki on migrating from ete3 to ete4. @lm-jkominek and @sarscov2 thanks and hope this helps for the moment!

The GitHub wiki was updated here. The Wiki page is kept up to date,

> "Jordi Burguet-Castell edited this page 4 days ago · 2 revisions"

The actual additional features (why ETE4 is needed) and how to access them are not listed here. However the changes in syntax are listed. Examples of changes ...

```
# ETE3
t = Tree(filename)
# ETE4
t = Tree(open(filename))
# ETE3
Tree(name='Hi', dist=1, support=1)
# ETE4
Tree({'name': 'Hi', 'dist': 1, 'support': 1})
# ETE3
Tree(nw, format=1)
# ETE4
Tree(nw, parser=1)

```

In summary there's different syntax. So we can start using ETE4

---
Tags: python, phylogenetics, phylogeny, ete3
---


# Question
Title: Looping over a list of cell models in Python
I have two data sets

1. Mutations in 10 genes
2. Gene expression in 18 cancer 3D models

*For every (mutation, gene) and gene expression pair we run an ANOVA test*

This is my data

```

> mut=read.delim("mut.tsv",row.names = 1)
> expr=as.data.frame(t(read.delim("expr.tsv",row.names = 1)))

dput(expr)

structure(list(GeneA_KO = c(0.0756383119639074, -0.162841965682723,
-0.110703100204517, -0.457154616239399, -0.505815193885352, 0.366072775033448,
0.288650458935083, -0.0287646498640309, -0.690732824311905, 0.157366122620995,
-1.29451546737684, 0.124167425460344, -0.168045313861751, 0.0186504943306465,
-0.495880827400706, -0.206668627339351, -0.148609490335447, 0.066844529089804
), GeneB_KO = c(-0.870157285953621, -0.35900872537626, -0.990406888128838,
-0.120656981579292, -1.41427631791049, -0.866421095836781, -0.816936112146099,
-0.894597236248507, -1.00084584207823, -0.850754488133456, -1.05314239761653,
-1.07580675428372, -1.24708324665688, -1.04098460460899, -0.232388675985086,
-1.48513703581865, -1.48373348870346, -1.76484058374354), class = "data.frame", row.names = c("ModelM",
"ModelO", "ModelS", "ModelA", "ModelI", "ModelC", "ModelT", "ModelH",
"ModelE", "ModelR", "ModelA2", "ModelP", "ModelE2", "ModelU",
"ModelT2", "ModelC2", "ModelS2", "ModelX"))

dput(mut)

structure(list(ModelM = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
0L), ModelO = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelS = c(0L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelA = c(0L, 0L, 0L, 1L,
0L, 0L, 0L, 0L, 0L, 0L), ModelI = c(0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 0L), ModelC = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L
), ModelT = c(0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L), ModelH = c(1L,
1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L), ModelE = c(1L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 1L), ModelR = c(1L, 0L, 0L, 1L, 0L, 0L, 1L,
0L, 1L, 0L), ModelA2 = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
0L), ModelP = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L), ModelE2 = c(1L,
0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L), ModelU = c(1L, 0L, 1L, 1L,
0L, 1L, 0L, 1L, 1L, 1L), ModelT2 = c(0L, 0L, 0L, 1L, 0L, 0L,
0L, 0L, 0L, 0L), ModelI2 = c(1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L,
0L, 0L), ModelC2 = c(1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L),
    ModelS2 = c(0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 1L)), class = "data.frame", row.names = c("Gene1_mut",
"Gene2_mut", "Gene3_mut", "Gene4_mut", "Gene5_mut", "Gene6_mut",
"Gene7_mut", "Gene8_mut", "Gene9_mut", "Gene10_mut"))

```

To test this, for the first 3D cancer model as instance, I have done like below

```
           GeneA_KO   GeneB_KO
ModelM   0.07563831 -0.8701573
ModelO  -0.16284197 -0.3590087
ModelS  -0.11070310 -0.9904069
ModelA  -0.45715462 -0.1206570
ModelI  -0.50581519 -1.4142763
ModelC   0.36607278 -0.8664211
ModelT   0.28865046 -0.8169361
ModelH  -0.02876465 -0.8945972
ModelE  -0.69073282 -1.0008458
ModelR   0.15736612 -0.8507545
ModelA2 -1.29451547 -1.0531424
ModelP   0.12416743 -1.0758068
ModelE2 -0.16804531 -1.2470832
ModelU   0.01865049 -1.0409846
ModelT2 -0.49588083 -0.2323887
ModelC2 -0.20666863 -1.4851370
ModelS2 -0.14860949 -1.4837335
ModelX   0.06684453 -1.7648406

```

```
           ModelM ModelO ModelS ModelA ModelI ModelC ModelT ModelH ModelE
Gene1_mut       0      0      0      0      0      0      0      1      1
Gene2_mut       0      0      0      0      0      0      0      1      0
Gene3_mut       0      0      0      0      0      0      0      1      0
Gene4_mut       1      1      1      1      1      1      1      1      0
Gene5_mut       0      0      0      0      0      0      0      0      1
Gene6_mut       0      0      0      0      0      0      1      1      0
Gene7_mut       0      0      0      0      0      0      0      1      0
Gene8_mut       0      0      0      0      0      0      1      1      0
Gene9_mut       0      0      0      0      0      0      1      0      0
Gene10_mut      0      0      0      0      0      0      0      0      1
           ModelR ModelA2 ModelP ModelE2 ModelU ModelT2 ModelI2 ModelC2 ModelS2
Gene1_mut       1       0      0       1      1       0       1       1       0
Gene2_mut       0       0      0       0      0       0       0       1       1
Gene3_mut       0       0      0       0      1       0       0       0       0
Gene4_mut       1       1      1       1      1       1       1       0       1
Gene5_mut       0       0      0       0      0       0       0       0       1
Gene6_mut       0       0      0       1      1       0       0       1       1
Gene7_mut       1       0      0       0      0       0       0       0       1
Gene8_mut       0       0      0       0      1       0       0       1       0
Gene9_mut       1       0      0       0      1       0       0       0       0
Gene10_mut      0       0      1       0      1       0       0       0       1

```

```
> my_data=cbind(mut[,1],expr[,1])
> colnames(my_data)=c("mut","expr")
> res.aov <- aov(expr ~ mut, data = as.data.frame(my_data))
> # Summary of the analysis
> summary(res.aov)
            Df Sum Sq Mean Sq F value Pr(>F)
mut          1  0.145  0.1448   0.192  0.673
Residuals    8  6.022  0.7527

```

Now I do a for loop for all of the 18 cancer cell lines

```
> for(i in 1:ncol(g))
+ {
+   my_data=as.data.frame(cbind(m[,1],g[,1]))
+   #tidy will summarise and return neat format
+   print(broom::tidy(aov(my_data[,2] ~ my_data[,1], data = my_data)))
+ }
# A tibble: 2 × 6
  term            df sumsq meansq statistic p.value
  <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
1 my_data[, 1]     1 0.145  0.145     0.192   0.673
2 Residuals        8 6.02   0.753    NA      NA
# A tibble: 2 × 6
  term            df sumsq meansq statistic p.value
  <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
1 my_data[, 1]     1 0.145  0.145     0.192   0.673
2 Residuals        8 6.02   0.753    NA      NA

```

I done this in R but of course I got when doing in Python

*Python*

```
import os
import pandas as pd
from statsmodels.formula.api import ols
import statsmodels.api as sm
os.getcwd()
os.chdir("/Users//Downloads/test")
os.getcwd()
m = pd.read_csv("mut.tsv", sep="\t", index_col=0,delim_whitespace=True)
m=pd.DataFrame(m)
m
g = pd.read_csv("expr.tsv", sep="\t",index_col=0,delim_whitespace=True).T
g=pd.DataFrame(g)
g
lst = list(g.columns)

lst
dfs = []
for Model in lst:
   dfs.append(pd.concat([m.iloc[:, 0].reset_index(drop=True), g.iloc[0,:].reset_index(drop=True)],axis=1))
   dfs.append(df.dropna())
   dfs.append(df.columns==['mut','expr'])
   model = ols('expr ~ C(mut)', data=df).fit()
   aov_table = sm.stats.anova_lm(model, typ=2)
   aov_table
   final = pd.concat(dfs, ignore_index=True)

```

```
    for Model in lst:
      Input In [165]
        for Model in lst:
                         ^
    IndentationError: expected an indented block

```

# Answer
Here is one potential approach:

Example data:

```
expr1 <- structure(list(GeneA_KO = c(0.0756383119639074, -0.162841965682723, -0.110703100204517, -0.457154616239399, -0.505815193885352, 0.366072775033448, 0.288650458935083, -0.0287646498640309, -0.690732824311905, 0.157366122620995, -1.29451546737684, 0.124167425460344, -0.168045313861751, 0.0186504943306465, -0.495880827400706, -0.206668627339351, -0.148609490335447, 0.066844529089804), GeneB_KO = c(-0.870157285953621, -0.35900872537626, -0.990406888128838, -0.120656981579292, -1.41427631791049, -0.866421095836781, -0.816936112146099, -0.894597236248507, -1.00084584207823, -0.850754488133456, -1.05314239761653, -1.07580675428372, -1.24708324665688, -1.04098460460899, -0.232388675985086, -1.48513703581865, -1.48373348870346, -1.76484058374354)), class = "data.frame", row.names = c("ModelM", "ModelO", "ModelS", "ModelA", "ModelI", "ModelC", "ModelT", "ModelH", "ModelE", "ModelR", "ModelA2", "ModelP", "ModelE2", "ModelU", "ModelT2", "ModelC2", "ModelS2", "ModelX"))
mut <- structure(list(ModelM = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelO = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelS = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelA = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelI = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelC = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelT = c(0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L), ModelH = c(1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L), ModelE = c(1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L), ModelR = c(1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L), ModelA2 = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelP = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L), ModelE2 = c(1L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L), ModelU = c(1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L), ModelT2 = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelI2 = c(1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelC2 = c(1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L),     ModelS2 = c(0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 1L)), class = "data.frame", row.names = c("Gene1_mut", "Gene2_mut", "Gene3_mut", "Gene4_mut", "Gene5_mut", "Gene6_mut", "Gene7_mut", "Gene8_mut", "Gene9_mut", "Gene10_mut"))

```

---

R code:

```
my_data=cbind(mut[,1],expr1[1:10,1])
colnames(my_data)=c("mut","expr1")
res.aov <- aov(expr1 ~ mut, data = as.data.frame(my_data))
summary(res.aov)
#>             Df Sum Sq Mean Sq F value Pr(>F)
#> mut          1 0.1364  0.1364   1.114  0.322
#> Residuals    8 0.9789  0.1224

for (j in 1:nrow(mut)) {
  for(i in 1:ncol(expr1)) {
    my_data=as.data.frame(cbind(mut[,j],expr1[1:10,i]))
    print(broom::tidy(aov(my_data[,2] ~ my_data[,1], data = my_data)))
  }
}
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df  sumsq meansq statistic p.value
#>   <chr>        <dbl>  <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.0612 0.0612     0.465   0.515
#> 2 Residuals        8 1.05   0.132     NA      NA
#> # A tibble: 2 × 6
#>   term            df  sumsq meansq statistic p.value
#>   <chr>        <dbl>  <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.0637 0.0637     0.479   0.508
#> 2 Residuals        8 1.06   0.133     NA      NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.246  0.246      2.26   0.171
#> 2 Residuals        8 0.869  0.109     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.313  0.313      3.07   0.118
#> 2 Residuals        8 0.815  0.102     NA     NA
#> # A tibble: 2 × 6
#>   term            df   sumsq  meansq statistic p.value
#>   <chr>        <dbl>   <dbl>   <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.00108 0.00108   0.00777   0.932
#> 2 Residuals        8 1.11    0.139    NA        NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.220  0.220      1.94   0.201
#> 2 Residuals        8 0.907  0.113     NA     NA
#> # A tibble: 2 × 6
#>   term            df  sumsq meansq statistic p.value
#>   <chr>        <dbl>  <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.0529 0.0529     0.398   0.546
#> 2 Residuals        8 1.06   0.133     NA      NA
#> # A tibble: 2 × 6
#>   term            df  sumsq meansq statistic p.value
#>   <chr>        <dbl>  <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.0901 0.0901     0.695   0.429
#> 2 Residuals        8 1.04   0.130     NA      NA

# write.csv(mut, "mut.csv", quote = FALSE, row.names = FALSE)
# write.csv(expr1[1:10,], "expr1.csv", quote = FALSE, row.names = FALSE)

```

Created on 2023-09-19 with reprex v2.0.2

---

Python code:

import os
import pandas as pd
from statsmodels.formula.api import ols
import statsmodels.api as sm
m = pd.read_csv("mut.csv", sep=",")
m=pd.DataFrame(m)
m.info()
g = pd.read_csv("expr1.csv", sep=",")
g=pd.DataFrame(g)

lst = list(g.columns)
m_lst = list(m.columns)

dfs = pd.DataFrame({})
for Mut in m_lst:
    for Model in lst:
        df = pd.concat([m[[Mut]].reset_index(drop=True), g[[Model]].reset_index(drop=True)],axis=1)
        df = df.rename(columns={Mut: "mut", Model: "g"})
        model = ols('g ~ C(mut)', data=df).fit()
        aov_table = sm.stats.anova_lm(model, typ=2)
        aov_table[['Model']] = Model
        aov_table[['Mut']] = Mut
        dfs = dfs.append(aov_table)

with pd.option_context('display.max_rows', None, 'display.max_columns', None):
    print(dfs)
            sum_sq   df         F    PR(>F)     Model      Mut
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelM
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelM
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelM
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelM
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelO
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelO
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelO
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelO
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelS
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelS
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelS
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelS
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelA
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelA
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelA
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelA
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelI
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelI
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelI
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelI
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelC
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelC
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelC
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelC
C(mut)    0.061205  1.0  0.464539  0.514753  GeneA_KO   ModelT
Residual  1.054037  8.0       NaN       NaN  GeneA_KO   ModelT
C(mut)    0.063734  1.0  0.479234  0.508361  GeneB_KO   ModelT
Residual  1.063933  8.0       NaN       NaN  GeneB_KO   ModelT
C(mut)    0.245964  1.0  2.263619  0.170857  GeneA_KO   ModelH
Residual  0.869278  8.0       NaN       NaN  GeneA_KO   ModelH
C(mut)    0.312937  1.0  3.072788  0.117706  GeneB_KO   ModelH
Residual  0.814730  8.0       NaN       NaN  GeneB_KO   ModelH
C(mut)    0.001082  1.0  0.007771  0.931920  GeneA_KO   ModelE
Residual  1.114159  8.0       NaN       NaN  GeneA_KO   ModelE
C(mut)    0.220171  1.0  1.940909  0.201060  GeneB_KO   ModelE
Residual  0.907496  8.0       NaN       NaN  GeneB_KO   ModelE
C(mut)    0.052891  1.0  0.398296  0.545575  GeneA_KO   ModelR
Residual  1.062350  8.0       NaN       NaN  GeneA_KO   ModelR
C(mut)    0.090105  1.0  0.694741  0.428739  GeneB_KO   ModelR
Residual  1.037562  8.0       NaN       NaN  GeneB_KO   ModelR
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO  ModelA2
Residual  0.978877  8.0       NaN       NaN  GeneA_KO  ModelA2
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO  ModelA2
Residual  0.586718  8.0       NaN       NaN  GeneB_KO  ModelA2
C(mut)    0.004637  1.0  0.033399  0.859537  GeneA_KO   ModelP
Residual  1.110605  8.0       NaN       NaN  GeneA_KO   ModelP
C(mut)    0.276724  1.0  2.601573  0.145422  GeneB_KO   ModelP
Residual  0.850943  8.0       NaN       NaN  GeneB_KO   ModelP
C(mut)    0.044310  1.0  0.331000  0.580888  GeneA_KO  ModelE2
Residual  1.070932  8.0       NaN       NaN  GeneA_KO  ModelE2
C(mut)    0.170278  1.0  1.422852  0.267110  GeneB_KO  ModelE2
Residual  0.957389  8.0       NaN       NaN  GeneB_KO  ModelE2
C(mut)    0.001687  1.0  0.012120  0.915049  GeneA_KO   ModelU
Residual  1.113555  8.0       NaN       NaN  GeneA_KO   ModelU
C(mut)    0.008679  1.0  0.062048  0.809563  GeneB_KO   ModelU
Residual  1.118988  8.0       NaN       NaN  GeneB_KO   ModelU
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO  ModelT2
Residual  0.978877  8.0       NaN       NaN  GeneA_KO  ModelT2
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO  ModelT2
Residual  0.586718  8.0       NaN       NaN  GeneB_KO  ModelT2
C(mut)    0.017610  1.0  0.128353  0.729416  GeneA_KO  ModelI2
Residual  1.097631  8.0       NaN       NaN  GeneA_KO  ModelI2
C(mut)    0.260821  1.0  2.407079  0.159387  GeneB_KO  ModelI2
Residual  0.866846  8.0       NaN       NaN  GeneB_KO  ModelI2
C(mut)    0.191207  1.0  1.655405  0.234206  GeneA_KO  ModelC2
Residual  0.924035  8.0       NaN       NaN  GeneA_KO  ModelC2
C(mut)    0.033474  1.0  0.244741  0.634101  GeneB_KO  ModelC2
Residual  1.094193  8.0       NaN       NaN  GeneB_KO  ModelC2
C(mut)    0.044621  1.0  0.333425  0.579533  GeneA_KO  ModelS2
Residual  1.070620  8.0       NaN       NaN  GeneA_KO  ModelS2
C(mut)    0.096956  1.0  0.752533  0.410940  GeneB_KO  ModelS2
Residual  1.030711  8.0       NaN       NaN  GeneB_KO  ModelS2
```

> 4 votes

---
Tags: python
---


# Question
Title: Do we design primers for single stranded DNA with same criteria as we do for double stranded DNA?
I am trying to design primers for a single stranded DNA virus. But I am having a dilemma as whether we design the primers same way as double stranded DNA?
do we have to design forward and reverse primers both or just a single primer forward or reverse works?

# Answer
> 3 votes
If you've e.g. a parvovirus you design primers as normal, i.e. must be a pair of primers as for double stranded DNA. Albeit there's only one strand present after the first round of PCR there both strands will be present and without the reverse primer the PCR can't work. If you want evidence of this, you've had a SARS-CoV-2 diagnosis? Well thats a forward and reverse primer PCR starting with ssRNA converted to cDNA, i.e. single copy and then a pair of primers PCR the resulting locus.

---
Tags: primer, pcr
---


# Question
Title: Changing transparency of markers in Chimera X
I am using the software Chimera X to visualize the 3d composition of molecules. I am very new to Chimera and only use it for visualization. The files I'm using are cmm files (xml-formatted) basically looking like below:

```
...
<marker id="14" x="0.018" y="-0.180" z="0.215" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="15" x="-0.032" y="-0.191" z="0.241" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="16" x="-0.033" y="-0.117" z="0.246" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="17" x="0.020" y="-0.138" z="0.244" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="18" x="0.047" y="-0.083" z="0.280" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="19" x="0.073" y="-0.061" z="0.240" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="20" x="0.010" y="-0.068" z="0.271" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
...
<link id1="0" id2="1" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="1" id2="2" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="2" id2="3" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="3" id2="4" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="4" id2="5" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
...

```

i.e. just a simple polymer consisting of alot of markers (spheres) and links between them (cylinders), and with colors for each one specified in rgb (see here). This is all well and good, but now I would like to set transparency for a subset of the markers and associated links. Preferably I would set this directly in the attributes in the cmm files, but I cannot find how to do this. The manual does mention that colours in general can be set with rgba where a means alpha level, but when trying to set a like seen in the xml block above, nothing happens (I tried values close to both 0 and 1; no difference)

The manual documents a function called 'transparency' which I tried to use with for example `transparency 50,a #0:10-1200` to make atoms 10-1200 50% transparent. But I just get the error: "Missing or invalid "percent" argument: Expected a number". If I only type `transparency 50` everything does turn transparent, however, I wouldn't really call it transparent since I cannot really see through them to see the other molecules on the other side, they just turn bleak.

Questions:

1. *(Would someone explain how I could simply change the actual transparency of a large set of markers and links)* (preferably by loading some file with specified marker id:s)? Also, is it possible to assign this directly in the cmm files, and if so, how?

UPDATE: The command `transparency #1:2492-end 90 'All'` worked to set the transparency of markers e.g. 2492-end. What I would still like to know is if it is possible to set this directly in the cmm-files.

# Answer
Ok, this is not an answer, but to better visualize the problem I recreated the question/answer on `Chimera` that has a different syntax for the `transparency` command see : https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/midas/transparency.html and https://www.cgl.ucsf.edu/chimerax/docs/user/commands/transparency.html .

Input file is `marker.cmm` , to be opened in `FILE MENU` of `VOLUME TRACER` from `VOLUME MENU` :

```
<marker_set name="marker set 1">
<marker id="14" x="0.018" y="-0.180" z="0.215" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="15" x="-0.032" y="-0.191" z="0.241" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="16" x="-0.033" y="-0.117" z="0.246" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="17" x="0.020" y="-0.138" z="0.244" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="18" x="0.047" y="-0.083" z="0.280" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="19" x="0.073" y="-0.061" z="0.240" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="20" x="0.010" y="-0.068" z="0.271" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<link id1="14" id2="15" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="15" id2="16" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="17" id2="18" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="18" id2="19" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="19" id2="20" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
</marker_set>

```

Image before command :

Image after `transparency 50,a,b #0:16-end`

(#0 being 0 model as the only loaded model on my `Chimera` instance) :

Better transparency visualizations :

Dont get why a='50' or a='0.5' in `marker`,`link` definition dont work with `chimera` too
> 0 votes

---
Tags: visualization, 3d-structure, software-usage
---