# Question
Title: Remove files after process is terminated
In my pipeline nextflow, I create a channel with channel.fromSRA and the channel contains lots of heavy files.fastq.gz. Then I have a first process to unzip files and transform them in files.fasta and a second process to count kmers. Is it possible to remove files.gz after first process and remove files.fasta after the second process. Can someone explain me how to do? They are too heavy to keep them in the work directory. Thanks.

In my pipeline I use kmer-ssr to count kmers present on my genome sample.

First I create a reads\_ch = Channel.fromSRA

In this channel there are 12 files (. fastq.gz) of 20 Gb each one.

Then the first process unzips them and transforms them in files.fasta.

Now in my work directory there are 12 more files of 80 Gb each one. So the work directory is now 1200 Gb.

The second process counts kmer present in reads and creates 12 files of 100 Gb each one. The work directory will be 2200 Gb.

The work directory is too heavy. I want to know how I can suppress the current file. fastq.gz after first process

# Answer
> 0 votes
Unfortunately, this remains an open issue with Nextflow: https://github.com/nextflow-io/nextflow/issues/452

However, writing large amounts of uncompressed data to disk is usually best avoided. If it can't be (e.g. maybe the file needs to be read multiple times), you'd need to merge the two processes so that the writing of the temporary file and the running of the tool is done in a single process. This would then let you remove the temporary input file as soon as it's no longer required.

Another way would be to pipe the uncompressed and FASTA formatted reads directly into your tool. Not all tools support this functionality obviously, but Kmer-SSR *can* indeed receive input from STDIN, which completely obviates the need to write any temporary files. Note that by default, Kmer-SSR writes also to STDOUT and this lets you write a compressed output file:

```
process kmer_ssr {

    cpus 4

    input:
    tuple val(sample), path(fastq)

    output:
    tuple val(sample), path("${sample}.tsv.gz")

    """
    zcat -f "${fastq}" \\
        | sed -n '1~4s/^@/>/p; 2~4p' \\
        | kmer-ssr -t ${task.cpus} \\
        | gzip \\
        > "${sample}.tsv.gz"
    """
}

```

Since the above uses a pipeline, make sure to add the `pipefail` to your command. One way is to simply apply it to all processes by adding the following to your `nextflow.config`:

```
process {

  shell = [ '/bin/bash', '-euo', 'pipefail' ]
}

```

From the `bash` man page:

> pipefail
>
> If set, the return value of a pipeline is the value of the last (rightmost) command to exit with a non-zero status, or zero if all commands in the pipeline exit successfully. This option is disabled by default.

---
Tags: fasta, fastq, nextflow, k-mer, workflow-management
---


title: How can I detect MCPyV in a FASTQ file?

text:
# Question
Title: How can I detect MCPyV in a FASTQ file?
I have received two NGS files from an NGS company, both are FASTQ files that correspond to reads from a tumor sample. I heavily suspect that MCPyV is present, and I am hoping to identify it.

What downloadable software can I use for viral-host integration detection? I have access to both Windows and Mac OS, but have no experience with Linux.

# Answer
You could start by trying out a metagenomics analysis program like kraken2. There are a few databases that combine bacteria, fungi, and viral searches.

One benefit of using kraken2 is that it has options to output the most likely taxa of matches at the sub-read level, which can help with identifying integrated or chimeric sequences.

I've written up the higher-level approach in more detail in this answer:

How do I measure the proportion of different bacteria in a sample from a high-throughput sequencer?

The sequence-level hypotheses about kmer origin can be found in the direct output of kraken2 (last column):

```
$ kraken2 --threads 10 --db /media/minion/minion_temp/db/kraken/metagenome --memory-mapping fastq_pass/barcode0?/PAK83314_pass_barcode0?_*.fastq.gz --report kraken2.report.txt > kraken2.out.txt
$ head kraken2.out.txt
C       a3e0b7c9-e6f3-40bf-900c-b1c59eb4977d    9606    208     0:70 623:3 0:53 9606:1 0:6 9606:32 131567:5 2759:4
C       8030d33f-3838-4ed3-bab6-5f3b79fd059e    575614  164     0:68 623:1 0:46 575614:2 0:13
C       5d87a648-799d-451f-b747-32a0fdf2517a    9606    190     0:108 9606:5 0:43
C       8aba667d-4a5b-4f8b-9d18-1b5f761b87b2    9606    267     0:124 9606:9 0:99 1:1
C       1b3501c2-6766-4ee8-b984-9e92470f37da    9606    335     0:70 28901:5 0:12 35885:1 0:23 9606:16 0:83 9606:17 0:6 9606:2 0:9 9606:5 0:11 9606:41
U       75c2625d-ab87-4614-8d0d-c9d12c7d1e19    0       242     0:208
U       902b0ec1-80d9-4e07-80c8-29a79191d931    0       171     0:137
C       56b7c8a8-4cfd-44ec-8922-62871b29d31c    9606    338     0:103 9606:5 0:22 9606:5 0:1 9606:13 0:55 9606:27 0:38 9606:3 0:32
C       68bf6718-2692-46e7-ba2d-f45952fd6138    9606    197     0:68 1:9 0:38 9606:19 0:29
C       81f0f4ae-3bb5-4b94-89ff-152ff3b6f5d7    9606    238     0:146 9606:46 0:12

```

I'm aware that you're interested in looking for integration of a *single* viral genome. This can be specifically looked at with kraken2 because it can handle custom-built indexes as well. However, the benefit of using a comprehensive index (e.g. the `MinusB` database from here) is that it will make it much easier to exclude the other numerous endogenous viruses that are already integrated into the genome, whereas using a MCPyV-only index may lead to false-positive matches to sequences that are from another virus.
> 0 votes

# Answer
FYI (non-OPs) it is very likely this sample has already been diagnosed and proven positive via a PCR test. I don't know the OP, but I am pretty familiar with oncovirus studies. To be clear, my guess is that PCR positive - doesn't mean cancer.

1. What I think is occurring (if I'm wrong hopefully the OP will say) - is the OP has a cancer patient(s) and has good evidence that the underlying cause is MCPyV: which is a very specific oncovirus.
2. What they want to do is identify viral integration in the patient genome as evidence of oncogenesis. To put this *general* field in perspective I think this issue is literally the rise and fall of many investigators and companies in cancer (again not focusing on a singular virus *per se* but generically on the oncoviruses). Keep in mind I have an explanatory bias towards viruses.
3. The molecular genetics of integration is likely to be complex.

In summary, the OP will particularly want to do (I suspect) is map the viral integration site. What is difficult is not diagnosis of the virus because at any time post-any sort of approximate assembly the genome can be locally `blasted` against the reference virus and its likely to have already been PCR'd. It's getting an assembly good enough to understanding the integration site - to demonstrate integration has occurred and thats is likely difficult (but I hope I'm wrong). Inverse PCR was used for this sort of thing - but thats a wetlab approach and PCR-based strategies with traditional Sanger sequencing can still be very useful here, but its laborious but was/is a traditional approach.

Thus, I suspect the OP is from a classical virological background venturing into NGS. I noted the OP used a human reference genome hg38 in a previous post.

If thats correct then ... the big issue in my opinion is obtaining a correct assembly of the putative virus-human genome **O** Steps would be,

* QC
* Assembly e.g. via `bwa-mem2` ... **O** **O**
* Construct a local `makeblastdb` database of the genome then ..
* `blastn` between the human genome and a viral query to identify the integration site
* Perform something like a `gatk` analysis to have a different perspective on the integration site via VCF
* hone in on the contig via `IGV` for optical analysis and study how the virus has integrated

**O** My eukaryotic genomics is outdated and the above approach is summarised in Genomics in the Cloud by Geraldine A. Van der Auwera, Brian D. O'Connor published through O'Reilly Media. I just focus on small genomes.

**O** **O** and **NOTE** I cannot remember if `bwa-mem2` uses a reference genome e.g. hg38 type thing (Auwera and O'Connor will describe it). If it does the assembly of the integration site, in my opinion, would be compromised and alternative assembly strategies might (probably would) be required. I don't think this will be trivial - I think it would need a skilled eukaryotic bioinformatician who could be told about the biology of viral integration. Keep in mind, eukaryotic genome assembly is not my thing.

However, riches await the investigator if this can be satisfactorily solved IMO.

---

To answer @gringers point, `Kraken2` would be useful for viral integration *per se* and yeah the results would be interesting because there's quite of them ... but in this case - I think - its an exact genetic sequence of a single species, thus I think direct blast-type blasting or `minimap2` mapping is a better strategy than `Kraken2`. Just my opinion.
> 0 votes

---
Tags: fastq, virus
---


# Question
Title: Correspondence of SARS-CoV-2 annotations (Nextstrain clades - Pango lineages)
Is there any annotation file providing correspondence between Nextclade and Pangolin variant nomenclature/annotations, to annotate some SARS-CoV-2 genomes from Gisaid with both these?

For the moment I could only find:

What I'm looking for is an annotation as follows:

```
lineage    clade
AY.43      21J
AY.4       21J
...        ...

```

The idea is to generate a table as follows:

```
genome                                 lineage    clade
hCoV-19/Germany/BW-RKI-I-195742/2021   AY.43      21J
...                                    ...        ...

```

The way I proceed now is the following:

* Parse the Pangolin's lineage designation file to find a target genome
* Look for the target genome in Gisaid (for quality check)
* Parse the metadata associated to the selected genome to get the lineage for verification (only Pangolin lineage is reported)

These steps allow me to build the first part of the table:

```
genome                                 lineage
hCoV-19/Germany/BW-RKI-I-195742/2021   AY.43
...                                    ...

```

Then:

* Paste the genome to Nextclade, wait for the analysis to be completed
* Report the calculated clade

With this information, I can then build the table:

```
genome                                 lineage    clade
hCoV-19/Germany/BW-RKI-I-195742/2021   AY.43      21J
...                                    ...        ...

```

The idea behind this is that both lineages and clades can be used for "tagging" a genome, depending on the aim (a given clade e.g. 21J includes many lineages e.g. AY.43 and AY.4 in the example above) and that I don't enjoy manual and error-prone procedures;) Thanks for your input.

---

The GISAID nomenclature via Nextstrain and overall tree is presented here

(Crossposted on Nextstrain's discord forum)

# Answer
A comment from a Nextclade's user satisfies my actual needs, I will post it here for future reference:

> Hi - one possible source is the metadata file for UCSC’s daily-updated tree of public SARS-CoV-2 sequences:

```
https://hgdownload.soe.ucsc.edu/goldenPath/wuhCor1/UShER_SARS-CoV-2/public-latest.metadata.tsv.gz

```

> It’s a big file – almost 3 million lines – because there are about 3M sequences in GenBank & COG-UK at this point and that number will probably grow. The last two columns of the file (10th & 11th) are the Nextstrain clade and Pango lineage assigned to each sequence based on its placement in the UCSC/UShER tree and the tree’s annotation of where each clade/lineage starts.

> This command will make a file relating each Nextstrain clade to at least one lineage (often to a lot of lineages, especially for Delta):

```
gunzip -c public-latest.metadata.tsv.gz | cut -f 10,11 | sort -u > cladeToLineages.tsv

```

> Some Pango lineages have no public sequences (only GISAID sequences, which can’t be shared at a public URL), so those lineages will be missing from that file. Some older lineages are not monophyletic on the UCSC/UShER tree, so those lineages are not annotated and will be missing. (In total, ~200 / ~1500 are missing.)
> 3 votes

# Answer
Remember not to omit the 'clarification' that in addition to Pango lineages, GISAID clades, Nextstrain clades there is also the WHO nomenclature. We could designate the `Bioinformatics StackExchange monophyly`, surely there is room for one more nomenclature?

Here's the WHO's response copied and uploaded from here

Just to mention NCBI has adopted the Pango lineage nomenclature available on all of its sequences available via this link

I thought that Omicron was Pango lineage B12, or perhaps B.1.2? Oh I see it has been reclassified following further sequence data.
> 0 votes

# Answer
This information is now easily available in JSON format here: https://github.com/corneliusroemer/pango-sequences/blob/main/data/pango-consensus-sequences\_summary.json. It's one big dictionary of dictionaries with Pango lineage names as outer keys and corresponding Nextstrain clades accessible via inner key:

```
{
 "A": {
  "lineage": "A",
  "nextstrainClade": "19B",
  "nucSubstitutions": [
   "C8782T",
   "T28144C"
  ],
  ...
 },
 "A.1": {
  "lineage": "A.1",
  "unaliased": "A.1",
  "nextstrainClade": "19B",
  ...
 },
 ...
}

```

It's straightforward to convert this to a TSV mapping using just minimal scripting.

Beware that not all Pango lineages pre-Omicron are contained therein. Small ones are missing currently.

Disclaimer: I created the repo.
> 0 votes

---
Tags: phylogenetics, sequence-annotation, covid-19, variants
---


# Question
Title: Why might parabricks not be writing output to the output file I provided?
This feels like a silly question, but I am testing parabricks bammetrics and nothing is being written to my output file. This is the command I am using:

```
docker run \
  --rm \
  --volume ~/workdir:/workdir \
  --volume ~/ouputdir:/outputdir \
  -w /workdir nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1 \
  pbrun \
  bammetrics \
  --ref /workdir/reference.fasta \
  --bam /workdir/cram.cram \
  --interval-file /workdir/regions.bed \
  --out-metrics-file /outputdir/metrics.txt

```

But metrics.txt doesn't have anything written in it when I open it after the command is run. Am I doing something wrong? I am very new to parabricks, so I assume I am missing something simple

I apologize if this is the wrong place to ask this question.

# Answer
> 3 votes
**EDIT Summary**

For future users who encounter the problem the solution is,

```
docker run --platform linux/amd64 --rm -it -u 0 --entrypoint /bin/sh  nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1

```

I am unhappy with the practice of limiting Docker access via blocking standard runtime: there are statements of intent on SE as a strategy to attempt to restrict Docker access. Parabricks isn't an isolated incident of deploying this setup within bioinformatics. However, I checked their forum and there is no statement of intent that it's deliberate. Parabrick's copyrights are here. There would be issues if there was ever an official statement to this extent. They would likely argue it was for other reasons to automate stuff.

---

*Original answer* ...

I don't know parabricks, but I know Docker. What you are attempting is to run Docker from the outside ... I would get on the inside to view the error log. Whats happening is the error is not being passed outside the container. Getting on the inside requires

```
docker run  --rm -it --name mydocker nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1

```

This *should* give you access to the container. If it doesn't post back and I'll forward a different command that will. It's impossible to block access to a Docker container albeit some people try.

For simplicity is worth trying

```
docker run  --rm -it --name mydocker clara-parabricks:4.0.1-1

```

Anyway, if you are using this approach you'll need,

```
docker cp *outside container* *insider container* nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1 and vice versa

```

... this is an amazing command and allows you to move data in and out of the container. You've got to get the data on the inside of the docker and `docker cp` does that very well. I would then look to running the whole thing from the inside of the container. Thats exactly how I run Dockers, then I exit and return when it's done. There might be a bit of clean via this approach, I'm not sure `--rm` might take care of it.

You'll also be able to look at the programs and copy it out onto hard drive. This last bit may be not advisable - the whole point of Docker is to protect your system from their algorithms ... but you can definitely do it and if you get it paired up with NextFlow thats cool.

I don't particularly want to pull the Docker container because they're often huge and I've no interest in the algorithm albeit I've spotted its at `nvcr.io/nvidia/clara/`

Docker, to be honest it's lazy by the developers. Docker is - techie packaging. So it's good for reviewers and behind web firewalls. It's useful but like you've just experienced if it flops you're stuck. The one positive thing about Docker (apart from DevOps) is it's generally tolerated to package someone else's algorithm as part of a pipeline and let a third party 'pull it'. This is because they'll be a Dockerfile which can reconstruct the whole thing without copyright infringement. So if that first developer objects ... the Docker developer would just issue the Dockerfile. `conda` and stuff max out on that concept.

---

There's a weird error as the OP described below.

**Right** Docker uses by default a tag `latest` that does not work for parabricks. The tag must be `4.0.1-1`. Could you retry the above code please?

---

**Update** For good measure try the following:

```
docker run  --rm -it --name mystuff  nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1 --entrypoint /bin/bash

```

Might not be bash, thus

```
docker run  --rm -it --name mystuff  nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1 --entrypoint /bin/sh

```

They shouldn't throw an error, but you might still find you are outside the container after running it. Lets see.

---

**Update**

The error looks okay, because a runtime error is exactly what would be expected if they are attempting to prevent access. Therefore please try:

```
docker run --platform linux/amd64 --rm -it -u 0 --entrypoint /bin/sh  nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1

```

This, or a variant of it, should work.

---
Tags: wgs, docker
---


# Question
Title: Funcotator reference file error in GATK4
I am using a function called Funcotator (here) in GATK4. The code I am trying to run is:

```
gatk Funcotator \
--variant chr21.vcf \
--reference hg38.fa.fai \
--ref-version hg38 \
--data-sources-path funcotator_dataSources.v1.2.20180329 \
--output chr21.funcotated.maf \
--output-file-format MAF

```

However, it gives me an error message like this:

```
java.lang.IllegalArgumentException: File is not a supported reference file type: /Users/data/Bam/hg38.fa.fai
at htsjdk.samtools.reference.ReferenceSequenceFileFactory.lambda$getFastaExtension$0(ReferenceSequenceFileFactory.java:253)
at java.util.Optional.orElseGet(Optional.java:267)
at htsjdk.samtools.reference.ReferenceSequenceFileFactory.getFastaExtension(ReferenceSequenceFileFactory.java:253)
at htsjdk.samtools.reference.ReferenceSequenceFileFactory.getDefaultDictionaryForReferenceSequence(ReferenceSequenceFileFactory.java:223)
at org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile.checkFastaPath(CachingIndexedFastaSequenceFile.java:180)
at org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile.<init>(CachingIndexedFastaSequenceFile.java:143)
at org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile.<init>(CachingIndexedFastaSequenceFile.java:125)
at org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile.<init>(CachingIndexedFastaSequenceFile.java:110)
at org.broadinstitute.hellbender.engine.ReferenceFileSource.<init>(ReferenceFileSource.java:35)
at org.broadinstitute.hellbender.engine.ReferenceDataSource.of(ReferenceDataSource.java:27)
at org.broadinstitute.hellbender.engine.GATKTool.initializeReference(GATKTool.java:439)
at org.broadinstitute.hellbender.engine.GATKTool.onStartup(GATKTool.java:722)
at org.broadinstitute.hellbender.engine.VariantWalker.onStartup(VariantWalker.java:45)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:138)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:192)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:211)
at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:160)
at org.broadinstitute.hellbender.Main.mainEntry(Main.java:203)
at org.broadinstitute.hellbender.Main.main(Main.java:289)

```

Does anyone know how to solve this problem? I tried to use the reference file downloaded from the GATK cloud, but it still doesn't work and has the same error message. Also, I wonder if the reference files must be the version provided by the GATK team? The one I am using now is not from the GATK team, and as they force me to provide `--data-sources-path`, I just copied the one they used on their post. Could that cause any problem in the further analysis? Thanks a lot ~~~

Here is an update: I first tried `hg38.fa` as my reference and got this error:

```
A USER ERROR has occurred: Fasta dict file file:///Users/data/Bam/hg38.dict for reference file:///Users/data/Bam/hg38.fa does not exist. Please see http://gatkforums.broadinstitute.org/discussion/1601/how-can-i-prepare-a-fasta-file-to-use-as-reference for help creating it.

```

So I generated the `hg38.fa.fai`, then I got the error message shown above. Then, I tried the reference file downloaded from the GATK cloud, `resources-broad-hg38-v0-Homo_sapiens_assembly38.fasta.fai`, but the same error appeared again.

Another update: I tried to use the `hg38.dict` as the reference and used the same code:

```
    gatk Funcotator \
    --variant chr21.vcf \
    --reference hg38.dict \
    --ref-version hg38 \
    --data-sources-path funcotator_dataSources.v1.2.20180329 \
    --output chr21.funcotated.maf \
    --output-file-format MAF

```

However, I got the same error message again:

```
java.lang.IllegalArgumentException: File is not a supported reference file type: /Users/data/hg38.dict

```

# Answer
You need to pass the fasta file instead of fa.fai or dict file as reference

```
    gatk Funcotator \
    --variant chr21.vcf \
    --reference hg38.fa \
    --ref-version hg38 \
    --data-sources-path funcotator_dataSources.v1.2.20180329 \
    --output chr21.funcotated.maf \
    --output-file-format MAF

```

The directory where you have hg38.fa must contain also the hg38.fa.fai and hg38.dict

Hope this helps

Zhan
> 1 votes

---
Tags: gatk
---


# Question
Title: Can we test the accuracy of Phred scores shown in FASTQ files
Recently I have ran a human WGS on the BGI DNBSEQ system, and their FASTQ quality scores seem to be quite impressive, where the Phred scores barely deteriorate along the read length when checked on FastQC.

Since the DNBSEQ technology is relatively new to me, I wonder if there's any bioinformatic approach to test whether the Phred scores listed in these FASTQ files are authentic or not? I'm kind of worried that their sequencer would just falsely assign high Phred scores to all bases.

# Answer
Yes, reported accuracy can be tested and compared with actual accuracy by mapping reads from a reference sample to a known-perfect reference sequence. It's usually done in aggregate, at a read level rather than at single-base level.

I like using Heng Li's gap-compressed identity statistic for this:

Alternatively, fish out a few reads at random, and use BLAST, which reports an identity score.

If the sample is not a reference sample with a precisely known sequence, then this approach will not work, because introduced variation will confound the comparison.

While it might be possible to create a new reference sequence from the reads, this introduces the possibility of systematic error in the sequencer that would be hidden by the new reference sequence, and when you're looking at error rates in the realm of q30 and above (as in your plot), those small differences are going to matter a lot.
> 3 votes

---
Tags: fastq, wgs, fastqc, phred-scores
---


# Question
Title: Should I compensate for two clusters of samples in microarray by hierarchical clustering and MDS?
I seem to have two groups of samples based on hierarchical clustering as well as MDS plot analysis using the sample label from the tree cut of the hierarchial clustering. I am analyzing microarray data. Should I adjust for the two different groups of samples in the model? For context, my covariate of interest is a continuous measure (air pollutant measure). This is a human cohort and not a case-control model and hence I was wondering if I should adjust for the different groups of samples. \[

MDS plot is here for the same samples:

# Answer
Very interesting, if I understand correctly you are asking whether a straight hierarchical clustering should account for a continuous variable because it is resulting in discrete groups (clusters) which seems contradictory for a continuous variable?

If that is the question the answer is strictly no. Do not use weighting for *this particular* analysis (other regression style analysis yes - but not this one). Thus 'weights' are for regression analysis not for multivariate statistics. This clustering analysis is about assessing the correlation and not trying to impose it - thats regression.

The 'weighting' that can be used is different for a tree based analysis - it accounts for rate variation within the expression data - you are assuming rate homogeneity. That doesn't always work - but it's complicated and your data looks good as is.

Moreover, I would suggest your have 4 clusters. What you want to do is map the continuous trait across all branches 'tips' of your clustering. This will break up your continuous variable (hopefully) and looking at the MDS very much into discrete groups. That is the result, in fact it may not be 2 clusters but 4 clusters - map all the 'tips' against the continuous variable. There're statistically randomisation tree methods to assess the robustness of the observation - you can also use bootstrapping.

The question is whether these discrete groups are in fact the best description of the data. The MDS suggests that the continuous variable "air pollution" is in continuous within gene expression - using a combination of clustering and MDS you highlight this facet of the data. There alternative analyses namely regression or GLM to assess the relationship between the continuous variable and your results - the MDS is suggestive of that. However, it's not a 'and', 'or', 'right', 'wrong' relationship a continuous variable can be broken up into discrete groups: if there is a 1:1 relationship between the gene expression and the continuous variable "air pollution" (MDS analysis) - that is a better description, if that relationship doesn't hold breaking into discrete clusters is the best description *thus far*.

Essentially this is statistical modelling and the question is what is the best model to fit all the data. I'm personally not a fan of hierarchical clustering as a tree description of the data - but in your data it looks good.

At the heart of the question you are asking whether there are boundary effects in the data as described by the hierarchical clustering? Thus is there a certain level of pollution (or whatever) causing one gene expression boundary and then another and so on ... Perfectly possible - but central to that hypothesis is whether clustering is the best description of the data rather than a regression analysis, i.e. is it continuous or is it discrete.

**Discrete vs. continuous** If it is discrete then k-means clustering would be a really, really good follow-up analysis to confirm the discrete "boundary" effects within the data. GLM is an analysis that might trap the data if it's continuous, there a lot of other methods in regression analysis. Basically if you map each "pollution" exposure against the MDS and if it's a perfect fit - trapping the result in a regression model outside eigenvalues would be the traditional thing to do. However, these days "unsupervised learning" (of which MDS is one method) is trendy so it may not be needed.

**k-means clustering** I would certainly consider k-means clustering in any case - this is because hierarchical clustering makes a rate assumption which might smash up the discrete groups. K-means can capture that without getting into complicated tree analysis. Basically, k-means you'll set 2 groups - assess the fit of the continuous variable, try 3 and then 4, how does the air-pollution variable fit? The data will also work well with supervised machine learning. However, it is genuinely continuous then regression analysis is best approach. With k-means you can also do parameterisation, i.e. whats the optimal number of groups represented by the data: if it comes back with '4' - which neatly segregate "air pollution" thats really good.

It's an exciting result. The way you write the result in the context of the conclusion is really important - i.e. the continuous verse discrete best fit analysis - in case you encounter a statistical modeller in peer-review.

My answer is a bit waffly, but is a representation of statistical modelling in this particular data, because two approaches may both be appropriate, but one is a better description than the other.

---

**Note**, if I've misunderstood the question and some sort of standardisation is needed between experiments - thats a different answer. Thus if you are saying that the groups represent different experiments, e.g. performed on different days or different arrays and there's risk of experimental artefact: thats about standardisation and the answer would be absolutely. Standardisation is essential to rule out experimental bias - particularly in microarray data.
> 1 votes

---
Tags: phylogenetics, phylogeny, microarray
---


# Question
Title: How can I identify the virulence factors in my FASTQ files generated from nanopore sequencing bacteria DNA?
I have performed genomic sequencing using nanopore technology and upon analysing using BLAST in the command line I have identified the presence of E.coli. I was wondering if I can add some code to identify the known virulence genes of E.coli which are on the NCBI BLAST database? Thank you.

# Answer
you could try to submit your sequences to VFDB: a reference database for bacterial virulence factors ; paper here : reference paper.

use the VFanalyzer tool

> VFanalyzer provides an automatic analysis pipeline for a systematic screen of known/potential VFs in given complete/draft bacterial genomes. Users can perform the analysis by one of the publicly available bacterial genomes either from our pre-downloaded genome list (Option 1) or via any valid GenBank accession(s) (Option 2). Or alternatively, users can also analyse any private genome data (raw FASTA sequence, predicted proteins, or annotated GenBank format genome) using Option 3. For query with raw FASTA genome, predicted CDSs will be identified using GLIMMER3 prior to the run of VFanalyzer. Please note that, VFanalyzer only accepts query genomes from the genera of bacterial pathogens with full information in VFDB, and it's designed for complete or nearly complete draft genomes. So you are at your own risk if the query genome are derived from other bacteria or only partially sequenced.

> NOTE: Due to connetion issue with NCBI, VFanalyzer can not fetch gbk file using Edirect correctly. Please download the file manually and use option3 to submit your jobs, we will find the other way to fix the problem ASAP.

Think you need to convert your FASTQ to FASTA first.

Other resource is : VirulenceFinder 2.0 or PathogenFinder 1.1; paper ref : In Silico Genotyping of Escherichia coli Isolates for Extraintestinal Virulence Genes by Use of Whole-Genome Sequencing Data, PathogenFinder - Distinguishing Friend from Foe Using Bacterial Whole Genome Sequence Data
> 1 votes

---
Tags: genomics, blast, bacteria, command-line
---


# Question
Title: How do I install shapeit4 on Ubuntu/Mint?
I'm having trouble install shapeit4 4.2.2 on my new Mint 21 installation (Ubuntu 22.04 LTS). According to the official docs, the requirements are:

* HTSlib: A great C library for reading/writing high-throughput sequencing data.
* BOOST: A free peer-reviewed portable C++ source libraries. SHAPEIT4 uses two specific BOOST libraries: iostreams and program\_options.

I installed HTSlib with `make` and `sudo make install`, and libboost using its calls `./bootstrap` and `./b2` but I still get errors when trying to run the `make` file for shapeit4:

```
/usr/bin/ld: cannot find /home/user/Tools/htslib-1.11/libhts.a: No such file or directory
/usr/bin/ld: cannot find /usr/lib/x86_64-linux-gnu/libboost_iostreams.a: No such file or directory
/usr/bin/ld: cannot find /usr/lib/x86_64-linux-gnu/libboost_program_options.a: No such file or directory

```

# Answer
Manually editing makefiles is just... ugh. Save yourself the effort and use conda:

```
conda create -n shapeit
conda install -c bioconda shapeit4 -n shapeit

```
> 1 votes

# Answer
The documentation appears to be missing an edit one has to do. Open the `makefile` in the shapeit4 dir, and you will see something like (for version 4.2.2):

```
#COMPILER MODE C++11
CXX=g++ -std=c++11

#HTSLIB LIBRARY [SPECIFY YOUR OWN PATHS]
HTSLIB_INC=$(HOME)/Tools/htslib-1.11
HTSLIB_LIB=$(HOME)/Tools/htslib-1.11/libhts.a

#BOOST IOSTREAM & PROGRAM_OPTION LIBRARIES [SPECIFY YOUR OWN PATHS]
BOOST_INC=/usr/include
BOOST_LIB_IO=/usr/lib/x86_64-linux-gnu/libboost_iostreams.a
BOOST_LIB_PO=/usr/lib/x86_64-linux-gnu/libboost_program_options.a

```

So the paths appear to be hardcoded. Set these to the paths where you installed htslib (and compile it first).

For libboost, I found that its own installer doesn't place files in this expected location. So you can either change these paths to where you installed libboost, or you can use the apt manager version which will put the files in the right location by default. This is done with `sudo apt-get install libboost-all-dev` (found here).

After this, just run the `make` file again, and the output should looks like this:

```
user@computer:/dir/shapeit4-4.2.2$ make
make: Nothing to be done for 'all'.

```
> 0 votes

# Answer
Actually, precompiled SHAPEIT 4.2.2 is available in a Ubuntu repository.

```
sudo apt install shapeit4

```
> 0 votes

---
Tags: software-installation, phasing
---


# Question
Title: Default placeholder for missing fields doesn't work with block in xtract tool
I am using the xtract tool of the Entrez Direct (EDirect) package of NBCI. I have a list of Accession Numbers:

```
U47804, U47803, U47802, U47801, U47800, U47799, U47798, X92938

```

and try to get information about:

1. Name of organism (Org-ref\_taxname tag)
2. Product of gene (Prot-ref\_name\_E tag)
3. Public year (Date-std\_year tag)

I executed the following command:

```
for i in {U47804 U47803 U47802 U47801 U47800 U47799 U47798 X92938}
do
  efetch -db "nucleotide" -id "$i" -format full; done |
    xtract -pattern Seq-entry -element Org-ref_taxname \
      -block Prot-ref_name -def "-" -element Prot-ref_name_E \
      -block Date-std -position first -element Date-std_year

```

The result I got:

```
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  1995

```

Because the **X92938** doesn't have *\<Prot-ref\_name\>*, so I was expecting to replace the missing fields with **"-"** by the `-def` tag of the extract tool. It means my desired result is as follows:

```
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1  gp120   1996
Human immunodeficiency virus 1    -     1995

```

Can someone have experience with the NBCI tool may help me? Thanks a lot.

# Answer
> 2 votes
After trying myself in many ways, I have found that the method of creating a checkmark column when a desired column is recorded is more efficient and easier to understand.

Here is my bash script to solve my problem in a verbose way (I guess has others way shorter).

```
#!/usr/bin/bash

change_index(){
    local in="$(< /dev/stdin)" # Variable assignment from the stdin
    local col_num=$(head -n 1 <(echo "$in") |awk -F"\t" '{print -1+NF/2}')  # Number of fields

    for k in $(seq 0 $(($col_num-1))); do # Loop through number of field
        in=$(awk -v ia="$k" -F"\t"  '{sub("/[^"ia"]","/"ia"\t-\t/"((ia+1)),$(2*ia+1)) ;print $0}' <(echo "$in")) # Supplement missing fields
    done

    in=$(awk -v ia="$col_num" -F"\t"  '{sub("/[^"ia"]","/"ia"\t-",$(2*ia+1)) ;print $0}' <(echo "$in")) # Supplement missing the last fields
    in=$(sed 's/\/[[:digit:]]\+//g' <(echo "$in")) # Remove markers
    echo "$in"
}

for i in {U47804 U47803 U47802 U47801 U47800 U47799 U47798 X92938}; do \
    efetch -db "nucleotide" -id "$i" -format full; done |
    xtract -pattern Seq-entry \
        -lbl "/0" -element Org-ref_taxname \
        -block Prot-ref_name -lbl "/1" -element Prot-ref_name_E \
        -block Date-std -position first -lbl "/2" -element Date-std_year |
    change_index| column -t

```

## Explanation

The main idea is that if a field is recorded there will be a marker field lay before the recorded field. This action is done by `-lbl` tag. And we have the following draft table:

```
/0  Human  immunodeficiency  virus  1  /1  gp120  /2  1996
/0  Human  immunodeficiency  virus  1  /1  gp120  /2  1996
/0  Human  immunodeficiency  virus  1  /1  gp120  /2  1996
/0  Human  immunodeficiency  virus  1  /1  gp120  /2  1996
/0  Human  immunodeficiency  virus  1  /1  gp120  /2  1996
/0  Human  immunodeficiency  virus  1  /1  gp120  /2  1996
/0  Human  immunodeficiency  virus  1  /1  gp120  /2  1996
/0  Human  immunodeficiency  virus  1  /2  1995

```

After having a draft table, and cleaning the table by `change_index()` function I wrote.

---
Tags: bash, edirect
---


# Question
Title: Python list comprehension to calculate the set of categories for each row
My data = `data`

```
Model   A1BG    A1CF    A2M A2ML1   A3GALT2 A4GALT  foldchange
IM00499 N   L   L   G   N   N   0.285744
IM00499.1   L   N   D   G   D   D   0.005046
IM00980 N   N   N   N   N   N   -0.068949
IM01529 N   N   L   N   N   N   0.065914
IM01624 A   N   L   N   N   N   0.302100
IM01240 N   D   N   D   N   G   -0.428324
IM00807 N   L   A   N   D   N   0.257222
IM00999 N   N   N   N   N   D   0.056268
IM00983 N   N   N   N   N   N   0.107456
IM01076 N   N   D   A   N   N   0.218301
IM00803 N   D   N   N   N   D   -0.479368
IM00125 N   N   A   N   N   A   0.270428
IM00998 D   N   N   N   G   N   0.142748
IM00997 N   L   G   G   N   G   -0.154540
IM00996 N   L   N   N   N   L   -1.302023

```

The above data frame is divided into 5 categories = A, N, L, G and D

`feature = ['A1BG', 'A1CF', 'A2M', 'A2ML1', 'A3GALT2', 'A4GALT']`

Feature\_level\_partitions are arrays of combination of classification values

`feature_level_partitions = ['A', 'N' , 'L', 'G', 'D'], ['D'], ['L', 'G']`

I was using the below `list comprehension` to calculate `Hedge's G`.

```
        hedges_g = np.array(
            [
                Hedges_g(data[feature].isin(feature_level_partition), y)
                for feature_level_partition in feature_level_partitions
            ],
            dtype=np.float64,
        )

```

`Hedges_g` is defined as follows:

```
def Hedges_g(x, y):
    difference_of_means = abs(np.mean(y[x])[0] - np.mean(y[~x])[0])
    s_pooled_weighted = math.sqrt((((n1 - 1) * s1 ** 2) + ((n2 - 1) * s2 ** 2)) / (n1 + n2 - 2))
    return difference_of_means / s_pooled_weighted

```

I am interested in knowing how many models there are within each category. So for example if I divide up the data into \["D","L"\] on one side, and the rest on the other, I want to know how many models have either "D" OR "L" and how many models have any of the other categories and then don't calculate `Hegde's G` if the number on either side is less than 2.

To answer the above query I tried the following code

```
[data[feature].isin(feature_level_partition).sum()
    for feature_level_partition in feature_level_partitions]

```

The above code gave me the number of models with the `feature_level_partitions` However I also want to know the number of models without the `feature_level_partitions`

I also want to create a python function which calculates both above numbers in one line and fit this code as a `if statement` in my above `list comprehension` function

# Answer
> 1 votes
What I'm not clear about is the definition of `y` \- thats crucial to understanding what the code does. Just to note, although you've called the dataframe `data`, I've called it `df` (traditional shorthand). In terms of speed there will not be much in it between vectorised `pandas` and `numpy`. Explicitly using `iterrows` or else `apply` is when it slows down a lot.

Possibilities of the definition of `y`.

1.

```
y = df.drop(['Model', 'foldchange'], axis=1)

```

I think the above is more likely, or ...

2.

```
y = df['fold change']

```

What I think is occurring is `np.mean` is performing a boolean operation on the feature columns. What I didn't understand though was this ...

```
abs(np.mean(y[x])[0] - np.mean(y[~x])[0])

```

.. this is because for example model 'IM00980' albeit the mean is 0 this will always be 1 under the above equation. This is because `y[~x]` will be 1 for all values and thus abs(0 - 1) is 1. Equally if all columns in a row were 1 the output would also be 1 because abs(1-0) is 1. Anyway, what I would do is explicitly perform the operation within `pandas` rather than coercing `np` into a boolean operation. This is

```
import pandas as pd

df = pd.read_csv(r'/pathtofile/megha.csv', delimiter="\t")

feature = ['A1BG', 'A1CF', 'A2M', 'A2ML1', 'A3GALT2', 'A4GALT']
df=df.drop(['foldchange'],axis=1)
feature_level_partitions = ['L', 'G']
filtered = df[feature].isin(feature_level_partitions).eq(True).mul(1)
df2 = df['Model'].to_frame().join(filtered)
df2['mean'] = df2[feature].mean(axis=1)
print(df2)

```

Here ~x would simply `1 - mean` as a new pandas column, i.e. `df2['~mean'] = df2['mean'] - 1`

The output of that code is:

| - | Model | A1BG | A1CF | A2M | A2ML1 | A3GALT2 | A4GALT | mean. |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 0 | IM00499 | 0 | 1. | 1. | 1. | 0. | 0. | 0.500000 |
| 1 | IM00499.1 | 1 | 0 | 0 | 1 | 0 | 0 | 0.333333 |
| 2 | IM00980 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000000 |
| 3 | IM01529 | 0 | 0 | 1 | 0 | 0 | 0 | 0.166667 |
| 4 | IM01624 | 0 | 0 | 1 | 0 | 0 | 0 | 0.166667 |
| 5 | IM01240 | 0 | 0 | 0 | 0 | 0 | 1 | 0.166667 |
| 6 | IM00807 | 0 | 1 | 0 | 0 | 0 | 0 | 0.166667 |
| 7 | IM00999 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000000 |
| 8 | IM00983 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000000 |
| 9 | IM01076 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000000 |
| 10 | IM00803 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000000 |
| 11 | IM00125 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000000 |
| 12 | IM00998 | 0 | 0 | 0 | 0 | 1 | 0 | 0.166667 |
| 13 | IM00997 | 0 | 1 | 1 | 1 | 0 | 1 | 0.666667 |
| 14 | IM00996 | 0 | 1 | 0 | 0 | 0 | 1 | 0.333333 |

---

If this is correct it provides the basis to move forward for all other operations. Basically, what I'm saying is I think is good to break open the `numpy` into explicit `pandas` operations, what is happening is there's a mix of `pandas` and `numpy` and I would recommend going down one route or the other.

---
Tags: python, pandas
---


# Question
Title: How do I build the "Standard Kraken 2 Database"?
I am at this point in the GitHub tutorial: https://github.com/DerrickWood/kraken2/blob/master/docs/MANUAL.markdown#standard-kraken-2-database

It says to create the standard Kraken database, I use the following command:

```
kraken2-build --standard --db $DBNAME

```

It then says to change `$DBNAME` with the preferred database name/location. Is this something that I need to download on my own? Where can I find it, and how do I do it?

To provide context, I want to use Kraken2 to identify whether or not there is the presence of MCPyV in a tumor sample.

# Answer
> 2 votes
\\$DBNAME is a folder location on your computer; you specify it yourself, for example:

```
kraken2-build --standard --db /media/ip80/externalDrive/kraken_indexes/betterThanGenbank

```

---
Tags: ngs, linux, metagenome
---


# Question
Title: How do I build the MinusB database for Kraken2? (Taxonomy issues)
I am attempting to build my own custom database for Kraken2. I have two questions:

1. If I have the MCPyV genome in a file called MCPyV.fasta, how do I build a database with just this?
2. How do I build the MinusB database? (https://benlangmead.github.io/aws-indexes/k2)

For number one, I attempted to build this by using:

```
kraken2-build --add-to-library MCPyV.fasta --db my_db

```

But was returned an error stating that there was no taxonomy information found. How do I make this MCPyV.fasta file amenable to Kraken2 while building the database?

For number 2, I again attempted to build the MinusB database using:

```
kraken2-build --build --threads 4 --db MinusB

```

and was returned an error regarding missing /taxonomy directory.

Does anyone know how to build these databases properly?

# Answer
I am not sure if I understand your first question correctly.

Kraken is meant to be for the taxonomic classification of metagenomic sequencing data. Therefore, if you only have one target genome of interest there is no need to use Kraken. You can just use any sequence aligner and check the mapping rates.

However, if you are trying to add a custom genome to an existing database then you need a taxonomy id for the sequence names in the genome and add it to your DB using the --add-to-library command. If this is what you are asking this tutorial would help you https://bioinformaticsworkbook.org/dataAnalysis/Metagenomics/Kraken.html#gsc.tab=0

Regarding your second question, I truly believe those ref indexes are ready to go. You can download the tar file, extract it in a folder called MinusB and release the Kraken using, e.g., `kraken2 --db MinusB yourreads.fq`

Hope this helps
> 1 votes

---
Tags: sequence-alignment, ngs, sequence-analysis, metagenome
---


# Question
Title: Phylogeny building in R from FASTA files:
Im working on building a phylogeny from scratch with downloaded FASTA sequences from GeneBank. I think Im doing alright up until the multi sequence alignment in the `msa` package where I get an error `Error in convertAlnRows(result$msa, type) : There is an invalid aln file!` As always thanks in advance for taking a look: Here is the reproducible example:

```
library (ape)
library (rentrez)
library (msa)

#Download Species Data
# B_terrestris
B_terrestris <- entrez_fetch(db = "nucleotide",
                                  id = "NC_045179.1",
                                  rettype = "fasta")

# B_hypocrita
B_hypocrita <- entrez_fetch(db = "nucleotide",
                             id = "NC_011923.1",
                             rettype = "fasta")

#Vespa
Vespa <- entrez_fetch(db = "nucleotide",
                      id = " MT137096.1",
                      rettype = "fasta")

seq<- c(B_terrestris, B_hypocrita, Vespa) # gotta clean it

#FASTA Cleaner Function
fasta_cleaner <- function(fasta_object, parse = TRUE)
{
  fasta_object <- sub("^(>)(.*?)(\\n)(.*)(\\n\\n)","\\4", fasta_object)
  fasta_object <- gsub("\n", "", fasta_object)

  if(parse == TRUE){
    fasta_object <- stringr::str_split(fasta_object,
                                       pattern = "",
                                       simplify = FALSE)
  }

  return(fasta_object[[1]])
}

#This should be ready to go!
clean.seq<-fasta_cleaner (seq, parse = FALSE)

str (clean.seq)

msa (clean.seq,  type="dna", method="Muscle", cluster="neighborjoining")

```

# Answer
The formal question I've been given is what aligner would I use? There is an issue between forming a data pipeline for `ape` and the latest and greatest. Its trade and the compromise would be `msaClustalOmega()`, but the rationale is complicated.

I strictly use `muscle5` or specifically `muscle -super5` option as a standalone. I am certain `R` `msa` has not introduced `muscle5` and I don't even know if it's published yet. It's on bioRxiv.org. The `super5` option must be specified and once done its amongst the most incredible algorithms *per se* notably in speed - no other aligner comes close - and personal benchmarks in accuracy. To give some idea, I was going to write an aligner as part of a package (its makes stuff simpler) and thought - whats the point with `muscle5`?

Its competitor was MAFFT and published benchmarks demonstrated it was the strongest, before `muscle5` was released.

The dilemma is simply that `ape` is cool and if you are wanting to have a seamless pipeline then greatest aligners don't fit with that. MAFFT isn't in `msa` and seriously doubt `muscle5` will be. HOWEVER Clustal Omega is in `msa` and that is the strongest (most recent) of the options it has. Omega does have published benchmarks and its neither quite as fast nor quite as accurate as MAFFT.

If you've got 50000 sequences - it's `muscle5` no contest. If you've 50 sequences then keep the pipeline together via `msaClustalOmega`.

---

The original question and code appeared to involve converting fasta inside `R` onto `msa`. Thats how subprocessing works outside `R` (or `Python`) but it's not the OO approach.

My original answer was in comments and remains an unofficial answer because usually I run the code I'm recommending and I don't fancy subprocessing muscle via `R`.

The error was the `msa` algorithm was attempting to shift the object into a format that will subprocess via muscle and it doesn't understand fasta format (it wants the raw object). Thus, you appear to be wanting to pass "official fasta" into msa. It doesn't need it, it simply needs a sequence object. It's imported as fasta but then it should be a sequence object (it's no longer fasta). That should go straight into msa. Thus `seq` should load straight into msa with "fasta" option removed. Finally could you remove msa and used `msaMuscle` then simply `msaMuscle(seq)`.
> 1 votes

---
Tags: r, phylogenetics, sequence-alignment, phylogeny, fasta
---


# Question
Title: Which candida albicans fasta and gff file should I use for alignment?
I downloaded the gff and fasta files for candida albicans from http://www.candidagenome.org. I want to use `hisat2` to align some fastq files against it.

There are many files in each download :

gff :

```
A22_ForcheSNPs.gff
A22_Historic_Assemblies.gff
A22_Jones_PMID_15123810_Polymorphisms.gff
A22_Jones_PMID_15123810_Polymorphisms.vcf
A22_Unannotated_transcripts_Bruno_et_al.gff
A22_Unannotated_transcripts_Sellam_et_al.gff
A22_Unannotated_transcripts_Tuch_et_al.gff
C_albicans_SC5314_A22_current_features.gff
C_albicans_SC5314_A22_current_features.gtf
C_albicans_SC5314_A22_current_features_with_chromosome_sequences.gff.gz
C_albicans_SC5314_A22_current_intergenic.gff
C_albicans_SC5314_haplotype_variations.gff
C_albicans_SC5314_version_A22-s07-m01-r179_features.gff
C_albicans_SC5314_version_A22-s07-m01-r179_features.gtf
C_albicans_SC5314_version_A22-s07-m01-r179_features_with_chromosome_sequences.gff.gz
C_albicans_SC5314_version_A22-s07-m01-r179_intergenic.gff

```

fasta :

```
C_albicans_SC5314_A22_current_chromosomes.fasta.gz
C_albicans_SC5314_A22_current_default_coding.fasta.gz
C_albicans_SC5314_A22_current_default_genomic.fasta.gz
C_albicans_SC5314_A22_current_default_protein.fasta.gz
C_albicans_SC5314_A22_current_not_feature.fasta.gz
C_albicans_SC5314_A22_current_orf_coding.fasta.gz
C_albicans_SC5314_A22_current_orf_genomic_1000.fasta.gz
C_albicans_SC5314_A22_current_orf_genomic.fasta.gz
C_albicans_SC5314_A22_current_orf_plus_intergenic.fasta.gz
C_albicans_SC5314_A22_current_orf_trans_all.fasta.gz
C_albicans_SC5314_A22_current_other_features_genomic_1000.fasta.gz
C_albicans_SC5314_A22_current_other_features_genomic.fasta.gz
C_albicans_SC5314_A22_current_other_features_no_introns.fasta.gz
C_albicans_SC5314_A22_current_other_features_plus_intergenic.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_chromosomes.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_default_coding.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_default_genomic.fasta
C_albicans_SC5314_version_A22-s07-m01-r179_default_protein.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_not_feature.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_orf_coding.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_orf_genomic_1000.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_orf_genomic.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_orf_plus_intergenic.fasta
C_albicans_SC5314_version_A22-s07-m01-r179_orf_trans_all.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_other_features_genomic_1000.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_other_features_genomic.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_other_features_no_introns.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_other_features_plus_intergenic.fasta.gz

```

I'm pretty sure I'd like to use the version `*A22-s07-m01-r179*` since it seems like it is a version that will is not transient like the `*current*` files.

QUESTION :

I am looking at bulk gene expression derived from an RNA-Seq experiment run on an Illumina sequencer. Which gff / fasta file should I use for alignment with `hisat`?

# Answer
First of all, you almost certainly want the `current` files. The convention in pretty much all sequence repositories is that `foo.current` is always the link pointing to the most recent and most up to date data. Those are what you want and not some "stable" version, you want the most accurate and that will usually be the `current`. In any case, my bet is that you will find that the `*A22-s07-m01-r179*` files are exactly the same as the corresponding `*current*` files since the `current` ones will just be symlinks.

Next, in order to align your data, you need the genome. The genome is this one:

C\_albicans\_SC5314\_A22\_current\_chromosomes.fasta.gz

That is the file that has the chromosomes in fasta format, so that's the one you want for alignment. Now, if you have done RNAseq and only care about annotated, known genes, you could speed things up a little by using one of the more targeted files (just scroll down on the page you linked to for an explanation of what each file has), but using the whole genome won't be wrong. And at 8.4M compressed, that's a tiny genome so I don't think it's worth trying to get much more speed out of what will already be a pretty quick process.

As for which GFF, that will again depend on what you want, but most likely you want either of these, depending on whether you want GFF or GTF:

```
C_albicans_SC5314_A22_current_features.gff
C_albicans_SC5314_A22_current_features.gtf

```

The main message is that you usually want the "current" files, just make a note of what version the "current" files are pointing to if you need to report that.
> 1 votes

# Answer
The refseq is *Candida albicans* SC5314.

I assume you are performing a fasta reference based assembly.

Its 8 chromosomes are NC\_032089.1 to NC\_032096.1 inclusively from chromosome 1 to chromosome 7 (NC\_032095.1) to chromosome R.

Its here.

Most of the files you downloaded are SC5314. So I dunno it depends what the project is re-annotation, reassembly?

The `vcf` is a very interesting file, particularly if you are sequencing your own variants of *Candida*.

Whats in all the files:

* C\_albicans\_SC5314\_A22\_current\_orf\_coding.fasta.gz proteins ignoring introns - useful
* C\_albicans\_SC5314\_A22\_current\_chromosomes.fasta.gz is going to be NCBI files I described.
* C\_albicans\_SC5314\_A22\_current\_default\_protein.fasta.gz orfs splicing out the introns

... list goes on really it depends on the project. I would recommend `current` only. `A22-s07-m01-r179` will likely be the previous assembly / annotation.

I guess the single file you want is C\_albicans\_SC5314\_A22\_current\_chromosomes.fasta.gz for referenced based assemblies (presumably of *Candida* isolates).
> 0 votes

---
Tags: rna-seq, fasta, gff
---


# Question
Title: computing the distance between ( XYZ coordinates)
I have a txt file (50k) proteins in it, which looks like this:

line 1:the protein code

line 2: protein length in amino acids

line 3: amino acid sequence

line 4: secondary structure

line 5: XYZ coordinates(which is numbers).

line 6:embty

the input looks like this :

```

4LGTD
247
M       S       E       K       L       Q       K       V       L       A       R       A       G       H       G       S       R   T
.       .       E       E       H       H       H       H       H       H       H       T       T       S       S       .       H   .
-34.3 -28.9 16.3        -33.0 -26.8 13.4        -29.4 -25.4 13.3        -27.9 -22.6 11.2        -26.9 -24.1 7.8 -23.3 -22.9 8.2 -23.0 -24.7 11.5     -24.1 -27.9 9.9 -21.8 -27.7 6.9 -18.8 -26.7 9.0 -19.5 -29.7 11.3        -19.6 -31.9 8.2 -16.2 -30.7 7.2 -16.7 -28.2 4.4      -15.0 -25.1 5.8 -15.4 -22.7 8.7 -18.8 -21.3 9.7 -17.7 -17.9 8.3 -16.3 -19.3 5.1 -19.5 -21.2 4.5 -21.5 -18.0 5.1 -19.7 -16.2 2.3      -20.8 -19.1 -0.0        -24.4 -19.2 1.2 -24.6 -15.4 0.6 -23.5 -16.0 -3.0        -26.4 -18.4 -3.6        -24.1 -21.4 -4.0    -25.9 -23.6 -1.5 -29.3 -25.2 -2.0        -31.4 -27.0 0.5 -34.0 -29.5 -0.8        -33.9 -27.8 -4.2        -34.1 -24.2 -3.0        -31.3 -21.7 -2.9     -30.7 -21.0 0.8 -30.2 -17.5 2.2 -27.9 -16.3 4.9 -29.4 -17.1 8.3 -31.3 -20.1 7.0 -31.7 -23.1 9.4 -32.1 -26.8 8.7 -32.6 -30.1 10.5     -30.2 -32.9 9.7 -32.0 -36.0 8.4 -31.6 -38.7 5.7 -34.0 -36.6 3.6 -32.0 -33.4 3.7 -30.4 -32.7 0.4 -27.7 -30.0 0.6 -25.8 -29.0 -2.5     -23.0 -26.6 -3.0        -22.5 -25.6 -6.6        -25.0 -28.3 -7.5        -22.9 -31.1 -5.9        -24.6 -33.1 -3.2    -23.0 -33.1 0.2  -23.1 -36.0 2.6 -23.7 -34.9 6.2 -23.8 -36.9 9.4 -27.2 -37.3 11.1        -28.5 -37.1 14.7        -19.2 -38.7 17.8    -17.7 -36.6 20.6 -14.0 -37.3 20.3        -11.2 -35.4 22.0        -7.4 -35.3 21.2 -4.7 -34.3 23.6 -1.9 -32.6 21.6 1.6 -31.6 22.7  2.2 -28.4 20.7       5.8 -27.5 20.1  7.1 -23.9 19.9  7.1 -22.6 16.3  3.6 -24.0 15.4  0.7 -21.6 14.7  -2.8 -22.2 16.0 -4.8 -21.4 12.8 -6.4 -23.3 9.9       -4.8 -20.9 7.5  -1.6 -18.8 8.1  -1.1 -16.0 5.5  2.7 -16.4 6.0   3.0 -20.1 5.2   4.2 -21.3 8.7   3.5 -24.7 10.3  0.1 -24.9 11.9       -1.1 -27.3 14.7 -3.6 -28.9 12.3 -1.0 -30.0 9.8  -0.7 -33.4 11.5 -4.0 -34.7 12.9 -6.1 -37.9 12.3 -8.5 -37.7 9.3  -12.2 -38.1 10.1     -15.3 -39.0 8.2 -17.5 -36.0 7.6 -16.2 -34.0 10.6        -13.2 -31.8 11.6        -10.9 -31.0 14.5        -11.7 -27.8 16.4     -8.8 -26.3 18.5 -9.8 -25.6 22.1 -8.4 -22.1 22.1 -5.0 -20.9 20.9 -1.5 -21.1 22.2 1.2 -18.6 21.2  3.8 -19.7 18.6  6.6 -20.0 21.1       4.5 -21.3 23.9  4.7 -25.1 24.5  2.3 -27.7 25.8  -1.3 -29.0 25.9 -3.8 -28.1 23.2 -7.0 -30.0 23.1 -8.9 -30.7 19.9 -12.5 -32.0 19.6     -14.9 -32.8 16.9        -17.6 -31.5 19.2        -17.9 -27.7 18.8        -20.2 -27.3 21.8        -17.6 -29.0 24.0        -15.0 -26.8 22.4     -16.8 -23.5 23.1        -17.7 -24.6 26.6        -14.1 -25.3 27.4        -12.9 -21.9 26.0        -15.6 -19.6 27.2    -14.8 -18.2 30.6 -18.4 -18.5 31.9        -17.5 -22.1 32.5        -14.3 -21.2 34.4        -12.2 -24.1 33.1        -8.7 -24.1 34.5 -5.7 -23.6 32.3      -2.5 -24.5 34.2 0.9 -23.2 33.2  4.4 -23.9 34.3  6.1 -20.6 34.0  9.8 -20.0 34.3  10.1 -16.4 35.4 13.0 -14.0 35.5 13.3 -11.2 38.1      12.7 -10.5 41.8 9.7 -12.2 43.3  8.4 -12.0 46.8  5.2 -13.0 48.6  3.7 -9.5 48.4   3.8 -9.8 44.6   2.2 -13.2 44.8  -0.6 -11.7 47.0      -1.0 -8.9 44.5  -1.3 -11.3 41.5 -3.9 -13.4 43.4 -6.1 -10.5 44.4 -6.1 -8.6 41.1  -3.9 -6.2 39.3  -4.5 -3.7 36.4  -3.6 -4.8 32.9       -3.8 -2.6 29.8  -7.0 -4.4 29.0  -8.5 -4.4 32.5  -8.2 -6.2 35.8  -6.8 -9.7 35.9  -5.7 -12.2 38.6 -4.2 -15.6 38.7 -6.6 -18.1 40.2      -3.8 -20.0 41.8 -0.0 -19.4 42.0 2.4 -21.9 43.4  6.2 -21.7 43.7  7.6 -24.7 41.8  11.4 -24.5 42.2 14.2 -22.0 42.2 17.4 -20.5 41.0      18.6 -17.0 41.7 16.7 -14.6 39.4 15.1 -17.6 37.8 12.2 -19.4 39.4 9.4 -21.8 38.3  5.8 -21.1 39.2  2.7 -23.1 38.4  -0.3 -20.8 37.8      -3.8 -21.1 36.7 -6.3 -18.9 34.9 -9.9 -19.4 34.1 -9.9 -17.1 31.1 -7.6 -16.7 28.1 -6.8 -13.1 27.2 -4.0 -12.5 24.7 -0.6 -13.1 26.3      -2.2 -12.9 29.7 0.3 -14.9 31.8  3.2 -12.8 30.5  1.3 -9.5 31.1   0.4 -10.5 34.6  4.0 -11.2 35.6  5.1 -8.0 33.8   2.5 -6.0 35.8        4.6 -6.8 38.9   8.0 -6.3 37.2   8.8 -9.9 36.2   9.7 -11.4 32.9  8.9 -14.8 31.5  11.6 -17.1 30.0 9.5 -20.2 29.5  5.8 -21.0 29.3       4.3 -24.5 29.1  0.6 -25.3 29.3  0.5 -28.5 31.4  -3.2 -28.9 32.0 -6.5 -27.9 30.5 -9.5 -28.7 32.7 -9.1 -32.2 34.0 -6.7 -33.3 31.2      -3.0 -32.9 32.1 -0.4 -33.0 29.4 -0.1 -36.8 28.8 2.9 -38.3 30.7  4.2 -40.6 28.0  3.5 -38.2 25.2  6.4 -36.4 23.5  6.5 -32.9 22.0       4.6 -32.8 18.6  2.5 -35.9 19.2  -1.1 -36.5 20.3 -3.4 -39.1 21.9 -7.1 -39.7 22.0 -9.2 -39.8 25.2 -11.5 -42.8 25.8    -15.2 -42.7 26.7 -14.8 -42.3 30.4        -12.3 -39.5 30.0        -14.2 -37.7 27.3        -17.3 -38.0 29.6        -15.4 -36.5 32.6    -14.3 -33.6 30.4 -17.9 -33.1 29.4        -19.2 -33.2 33.0        -16.7 -30.4 33.8        -18.7 -28.1 31.6        -22.2 -29.4 32.6    -22.6 -31.6 29.5 -24.2 -35.0 30.3        -22.3 -38.2 29.6        -22.5 -39.9 26.3        -22.8 -43.7 25.4

1XXIJ
334
M       R       W       Y       P       W       L       R       P       D       F       E       K       L       V       A       S   L
.       .       .       .       T       T       H       H       H       H       H       H       H       H       H       T       T   .
65.7 -16.7 128.5        64.5 -19.9 130.3        66.9 -22.7 131.5        66.0 -25.8 129.7        67.1 -29.4 130.9        69.2 -30.4 127.8     71.4 -27.4 128.6        72.9 -29.3 131.5        75.3 -32.0 130.1        76.9 -29.2 128.1        76.9 -26.7 131.0        78.6 -29.4 133.1     81.8 -29.7 131.1        81.8 -26.0 130.4        82.2 -25.2 134.1        84.3 -28.3 134.6        87.0 -26.5 132.5    86.5 -23.1 134.2 87.2 -24.8 137.5        90.3 -26.9 136.4        91.9 -24.1 134.3        91.2 -25.8 130.9        88.8 -22.9 130.1    88.9 -21.4 126.6 89.0 -17.6 126.7        86.9 -17.5 123.5        83.8 -19.6 123.0        81.6 -19.9 119.9        78.4 -21.8 120.5    76.4 -23.0 117.5 72.8 -23.9 117.9        69.4 -23.9 116.3        66.9 -21.3 117.2        65.1 -21.8 120.5        67.7 -24.4 121.7    68.4 -22.5 124.8 71.4 -20.9 123.2        70.2 -17.9 125.1        70.3 -19.4 128.5        73.6 -21.3 128.3        75.7 -18.2 127.8    74.1 -16.3 130.7 74.1 -19.4 133.0        77.8 -19.8 132.3        78.8 -16.2 133.0        77.0 -16.2 136.4        79.4 -18.9 137.5    82.4 -16.9 136.6 81.1 -13.6 138.0        80.6 -14.9 141.5        82.7 -13.6 144.4        81.9 -16.6 146.7        80.7 -19.6 144.5    79.2 -22.6 146.2 80.7 -25.3 143.9        78.7 -25.6 140.7        76.6 -22.5 141.0        76.6 -19.1 142.8        75.6 -18.0 146.3    73.1 -15.2 145.4 74.5 -12.7 148.0        78.0 -11.6 146.9        78.5 -8.3 145.0 78.7 -10.2 141.8        75.6 -12.3 142.3        73.7 -9.2 143.5      74.8 -6.8 140.8 74.2 -9.8 138.5 70.6 -10.2 139.7        70.4 -6.4 139.2 71.7 -7.1 135.7 74.0 -4.0 135.8 77.6 -5.1 136.4      81.1 -3.7 135.7 82.4 -7.0 134.2 79.5 -8.2 132.2 79.4 -6.7 128.7 76.5 -7.3 126.2 77.5 -6.8 122.6 74.1 -6.4 121.0 74.1 -5.1 117.4      73.4 -1.3 116.8 69.7 -0.7 117.7 67.9 -0.2 114.4 70.3 -1.8 112.0 70.7 -5.6 111.3 74.4 -5.6 111.8 77.1 -4.7 114.4 79.7 -1.9 114.3      83.5 -1.7 114.3 82.9 1.0 117.0  81.1 -1.0 119.7 83.9 -3.6 119.5 86.8 -1.2 119.9 85.1 0.5 122.9  84.7 -2.8 124.7 88.3 -3.9 123.7      89.8 -0.9 125.5 87.2 -1.3 128.2 88.8 -4.6 129.1 92.3 -3.2 128.8 91.9 -1.8 132.3 91.6 -3.7 135.7 88.3 -5.5 136.6 86.3 -3.2 138.9      85.3 -6.1 141.2 88.7 -7.8 141.0 88.2 -11.1 139.2        87.1 -12.3 135.8        85.6 -10.1 132.9        83.1 -11.6 130.4    81.8 -10.7 126.9 78.9 -12.0 124.9        77.8 -11.6 121.3        74.2 -12.0 120.1        75.2 -13.4 116.7        78.7 -13.1 115.2    77.5 -13.8 111.7 75.7 -10.4 111.7        78.8 -8.2 112.1 80.6 -5.9 109.6 83.7 -7.5 107.9 85.8 -4.8 109.5 84.2 -5.5 112.9 85.6 -9.1 112.6      89.2 -8.1 112.3 88.8 -5.4 114.9 87.8 -8.1 117.3 90.6 -10.6 116.4        93.6 -8.6 117.7 92.2 -8.7 121.2 91.9 -12.5 121.2    95.5 -12.5 120.1 96.5 -10.2 123.0        94.2 -11.8 125.7        93.9 -9.5 128.8 94.7 -11.9 131.7        91.7 -13.9 133.1        89.1 -12.4 130.8     86.3 -14.5 129.1        84.5 -14.4 125.7        81.3 -15.9 124.4        80.5 -15.9 120.7        77.3 -17.4 119.5    76.0 -18.2 116.0 73.5 -20.3 114.2        75.8 -20.6 111.2        79.5 -20.8 112.0        80.6 -21.2 108.4        79.5 -17.8 107.3    82.1 -16.6 109.9 85.7 -15.6 109.1        87.9 -18.7 109.1        90.7 -16.8 110.8        88.2 -15.8 113.5        86.5 -19.2 114.1    90.0 -20.6 114.6 90.7 -18.5 117.7        87.6 -19.7 119.5        86.5 -23.1 120.6        83.4 -24.2 118.7        80.9 -26.1 120.8    77.9 -27.9 119.3 74.6 -27.5 121.1        73.0 -30.7 119.6        69.2 -30.8 118.8        68.0 -33.4 121.2        66.8 -36.7 119.7    63.3 -36.2 118.3 61.6 -38.6 120.7        62.8 -36.6 123.6        61.9 -33.2 122.0        58.4 -34.3 121.2        57.6 -35.6 124.7    59.0 -32.4 126.2 56.8 -30.5 123.8        53.6 -32.5 124.5        53.8 -30.9 127.9        54.2 -27.3 126.9        51.2 -27.5 124.5    48.6 -30.1 123.9 48.5 -30.9 120.2        48.7 -34.1 118.1        52.2 -35.6 118.2        52.1 -35.5 114.4        51.8 -31.8 114.5    54.3 -31.6 117.3 56.9 -33.4 115.1        56.1 -31.1 112.2        56.9 -28.0 114.3        60.2 -29.7 115.5        61.2 -30.5 111.9    60.2 -27.0 110.8 62.8 -25.7 113.3        65.7 -27.9 112.4        65.1 -29.8 115.6        65.1 -26.7 117.9        62.5 -27.5 120.5    62.4 -23.9 121.6 61.6 -22.5 118.2        59.2 -25.4 117.7        57.4 -24.4 120.9        57.2 -20.8 119.7        55.7 -22.0 116.5    52.5 -22.6 118.4 52.3 -19.7 120.9        51.0 -17.2 118.3        48.2 -18.8 116.5        50.9 -19.8 114.0        50.0 -23.4 113.8    46.4 -22.5 113.1 47.8 -20.1 110.4        49.7 -23.1 109.1        46.6 -25.3 109.0        45.0 -22.6 106.8        47.8 -23.0 104.4    47.6 -26.8 104.0 43.8 -26.5 103.4        44.3 -24.0 100.7        47.1 -26.1 99.2 44.9 -29.2 99.2 42.3 -27.4 97.2 44.7 -26.0 94.7 46.9 -29.1 94.1      44.1 -31.1 92.5 42.2 -28.1 91.1 45.1 -26.2 89.4 43.8 -22.9 90.9 47.0 -22.0 92.7 46.0 -18.4 92.8 43.6 -19.4 95.6 46.6 -19.4 98.0      46.7 -15.7 97.7 44.6 -15.2 100.8        46.6 -17.8 102.7        50.0 -16.3 101.8        49.4 -12.6 101.4        50.1 -11.5 104.9     53.0 -9.8 106.8 54.0 -13.4 107.8        54.5 -14.0 104.1        58.1 -15.2 104.4        57.5 -17.6 107.2        54.6 -18.9 105.2     56.4 -19.2 101.9        59.0 -21.3 103.7        56.2 -23.6 105.0        55.0 -24.0 101.4        58.6 -24.8 100.4    59.0 -27.6 102.9 56.0 -29.3 101.4        57.5 -28.8 97.9 60.7 -30.3 99.0 59.0 -33.4 100.5        57.0 -34.1 97.3 60.3 -33.9 95.5 62.0 -36.3 97.9      59.2 -38.8 97.5 59.9 -39.0 93.7 63.3 -40.3 94.6 61.9 -42.6 97.2 63.3 -40.7 100.2        62.2 -42.2 103.5        62.3 -39.0 105.6     59.6 -36.3 105.3        59.3 -33.1 107.2        55.6 -32.7 106.8        54.4 -35.9 108.5        51.5 -33.8 109.6    49.8 -34.3 106.2 51.3 -37.2 104.2        48.2 -37.2 102.1        49.2 -33.8 100.8        52.8 -34.5 99.9 51.5 -37.6 98.1 48.9 -36.0 95.9      51.3 -33.1 95.3 54.0 -35.5 94.1 51.3 -37.3 92.0 50.2 -34.2 90.0 53.6 -32.6 89.7 56.6 -34.0 87.9 59.9 -33.8 89.7 60.9 -31.0 87.3      57.8 -28.9 87.8 58.3 -29.2 91.5 62.0 -28.0 91.2 60.8 -25.0 89.2 58.2 -24.0 91.7 60.5 -24.7 94.7 62.8 -22.0 93.4 60.4 -19.2 92.5      58.9 -19.8 95.9 62.4 -19.6 97.2 63.7 -16.5 95.4 60.4 -14.7 95.9 60.3 -15.3 99.6 63.7 -13.7 99.7 62.8 -10.5 97.9 60.1 -10.1 100.5     62.5 -10.4 103.3        65.4 -8.5 101.8 64.0 -5.6 99.5  62.8 -2.7 101.7 59.2 -2.6 100.6 57.9 -4.3 97.4  54.2 -5.0 96.9  53.5 -8.2 98.9       50.4 -8.5 96.9  51.8 -8.1 93.4  54.8 -10.2 93.8 52.9 -13.1 95.2 50.1 -12.8 92.6 52.9 -13.0 90.1 54.5 -16.0 91.9 51.1 -17.7 91.9      50.0 -17.1 88.3 53.5 -18.0 87.3 53.3 -21.3 89.1 50.3 -22.1 87.0 52.1 -21.6 83.7 54.8 -24.0 84.7 52.0 -26.5 85.2 50.8 -26.1 81.6      52.0 -28.7 79.0 54.8 -27.5 76.7 56.2 -24.8 78.9 60.0 -24.6 79.4 60.6 -24.0 83.1 63.6 -22.0 84.5 67.2 -23.0 85.4 69.4 -21.4 88.1      71.8 -18.6 87.1 75.6 -17.8 87.3

```

and so on (the same structure)

Please note in line 5; There is a space between x and y and y and z, but each XYZ group is tab-separated from the next.

Now I am unsure if the code down works and computing the distance between ( XYZ coordinates) correctly or not!

```
awk '
NR%6==5 { printf "%f\n", sqrt(($1-x)^2 + ($2-y)^2 + ($3-z)^2) }
       { x=$1; y=$2; z=$3 }' myfile.adataset

```

But I want the output to be like this, which is to print also the lines 1,2,3,4 , and the results:

```
4LGTD
247
M       S       E       K       L       Q       K       V       L
.       .       E       E       H       H       H       H       H
the result

1XXIJ
334
M       R       W       Y       P       W       L       R       P
        T       T       H       H       H       H       H
the result

and so on

```

So after that, I want other code to delete blocks of lines(1,2,3,4,5)if the results of XYZ coordinates are bigger than (4.5)? Any suggestions code for first computing the distance and then deleting if the line 5 result is bigger than 4.5 by awk to do this?

# Answer
Based on the example you've provided, your awk command is returning an incorrect result (XY problem), however one potential answer to your question is to print your desired outcome using two clauses:

`NR%6==0 (i.e. every 6 lines: 5 lines and 1 blank line) {do action}`

and

`NR%6!=0 (i.e. every other line) {do different action}`.

---

Example data:

```
cat test.txt

4LGTD
247
M       S       E       K       L       Q       K       V       L       A       R       A       G       H       G       S       R   T
.       .       E       E       H       H       H       H       H       H       H       T       T       S       S       .       H   .
-34.3 -28.9 16.3        -33.0 -26.8 13.4        -29.4 -25.4 13.3        -27.9 -22.6 11.2        -26.9 -24.1 7.8 -23.3 -22.9 8.2 -23.0 -24.7 11.5     -24.1 -27.9 9.9 -21.8 -27.7 6.9 -18.8 -26.7 9.0 -19.5 -29.7 11.3        -19.6 -31.9 8.2 -16.2 -30.7 7.2 -16.7 -28.2 4.4      -15.0 -25.1 5.8 -15.4 -22.7 8.7 -18.8 -21.3 9.7 -17.7 -17.9 8.3 -16.3 -19.3 5.1 -19.5 -21.2 4.5 -21.5 -18.0 5.1 -19.7 -16.2 2.3      -20.8 -19.1 -0.0        -24.4 -19.2 1.2 -24.6 -15.4 0.6 -23.5 -16.0 -3.0        -26.4 -18.4 -3.6        -24.1 -21.4 -4.0    -25.9 -23.6 -1.5 -29.3 -25.2 -2.0        -31.4 -27.0 0.5 -34.0 -29.5 -0.8        -33.9 -27.8 -4.2        -34.1 -24.2 -3.0        -31.3 -21.7 -2.9     -30.7 -21.0 0.8 -30.2 -17.5 2.2 -27.9 -16.3 4.9 -29.4 -17.1 8.3 -31.3 -20.1 7.0 -31.7 -23.1 9.4 -32.1 -26.8 8.7 -32.6 -30.1 10.5     -30.2 -32.9 9.7 -32.0 -36.0 8.4 -31.6 -38.7 5.7 -34.0 -36.6 3.6 -32.0 -33.4 3.7 -30.4 -32.7 0.4 -27.7 -30.0 0.6 -25.8 -29.0 -2.5     -23.0 -26.6 -3.0        -22.5 -25.6 -6.6        -25.0 -28.3 -7.5        -22.9 -31.1 -5.9        -24.6 -33.1 -3.2    -23.0 -33.1 0.2  -23.1 -36.0 2.6 -23.7 -34.9 6.2 -23.8 -36.9 9.4 -27.2 -37.3 11.1        -28.5 -37.1 14.7        -19.2 -38.7 17.8    -17.7 -36.6 20.6 -14.0 -37.3 20.3        -11.2 -35.4 22.0        -7.4 -35.3 21.2 -4.7 -34.3 23.6 -1.9 -32.6 21.6 1.6 -31.6 22.7  2.2 -28.4 20.7       5.8 -27.5 20.1  7.1 -23.9 19.9  7.1 -22.6 16.3  3.6 -24.0 15.4  0.7 -21.6 14.7  -2.8 -22.2 16.0 -4.8 -21.4 12.8 -6.4 -23.3 9.9       -4.8 -20.9 7.5  -1.6 -18.8 8.1  -1.1 -16.0 5.5  2.7 -16.4 6.0   3.0 -20.1 5.2   4.2 -21.3 8.7   3.5 -24.7 10.3  0.1 -24.9 11.9       -1.1 -27.3 14.7 -3.6 -28.9 12.3 -1.0 -30.0 9.8  -0.7 -33.4 11.5 -4.0 -34.7 12.9 -6.1 -37.9 12.3 -8.5 -37.7 9.3  -12.2 -38.1 10.1     -15.3 -39.0 8.2 -17.5 -36.0 7.6 -16.2 -34.0 10.6        -13.2 -31.8 11.6        -10.9 -31.0 14.5        -11.7 -27.8 16.4     -8.8 -26.3 18.5 -9.8 -25.6 22.1 -8.4 -22.1 22.1 -5.0 -20.9 20.9 -1.5 -21.1 22.2 1.2 -18.6 21.2  3.8 -19.7 18.6  6.6 -20.0 21.1       4.5 -21.3 23.9  4.7 -25.1 24.5  2.3 -27.7 25.8  -1.3 -29.0 25.9 -3.8 -28.1 23.2 -7.0 -30.0 23.1 -8.9 -30.7 19.9 -12.5 -32.0 19.6     -14.9 -32.8 16.9        -17.6 -31.5 19.2        -17.9 -27.7 18.8        -20.2 -27.3 21.8        -17.6 -29.0 24.0        -15.0 -26.8 22.4     -16.8 -23.5 23.1        -17.7 -24.6 26.6        -14.1 -25.3 27.4        -12.9 -21.9 26.0        -15.6 -19.6 27.2    -14.8 -18.2 30.6 -18.4 -18.5 31.9        -17.5 -22.1 32.5        -14.3 -21.2 34.4        -12.2 -24.1 33.1        -8.7 -24.1 34.5 -5.7 -23.6 32.3      -2.5 -24.5 34.2 0.9 -23.2 33.2  4.4 -23.9 34.3  6.1 -20.6 34.0  9.8 -20.0 34.3  10.1 -16.4 35.4 13.0 -14.0 35.5 13.3 -11.2 38.1      12.7 -10.5 41.8 9.7 -12.2 43.3  8.4 -12.0 46.8  5.2 -13.0 48.6  3.7 -9.5 48.4   3.8 -9.8 44.6   2.2 -13.2 44.8  -0.6 -11.7 47.0      -1.0 -8.9 44.5  -1.3 -11.3 41.5 -3.9 -13.4 43.4 -6.1 -10.5 44.4 -6.1 -8.6 41.1  -3.9 -6.2 39.3  -4.5 -3.7 36.4  -3.6 -4.8 32.9       -3.8 -2.6 29.8  -7.0 -4.4 29.0  -8.5 -4.4 32.5  -8.2 -6.2 35.8  -6.8 -9.7 35.9  -5.7 -12.2 38.6 -4.2 -15.6 38.7 -6.6 -18.1 40.2      -3.8 -20.0 41.8 -0.0 -19.4 42.0 2.4 -21.9 43.4  6.2 -21.7 43.7  7.6 -24.7 41.8  11.4 -24.5 42.2 14.2 -22.0 42.2 17.4 -20.5 41.0      18.6 -17.0 41.7 16.7 -14.6 39.4 15.1 -17.6 37.8 12.2 -19.4 39.4 9.4 -21.8 38.3  5.8 -21.1 39.2  2.7 -23.1 38.4  -0.3 -20.8 37.8      -3.8 -21.1 36.7 -6.3 -18.9 34.9 -9.9 -19.4 34.1 -9.9 -17.1 31.1 -7.6 -16.7 28.1 -6.8 -13.1 27.2 -4.0 -12.5 24.7 -0.6 -13.1 26.3      -2.2 -12.9 29.7 0.3 -14.9 31.8  3.2 -12.8 30.5  1.3 -9.5 31.1   0.4 -10.5 34.6  4.0 -11.2 35.6  5.1 -8.0 33.8   2.5 -6.0 35.8        4.6 -6.8 38.9   8.0 -6.3 37.2   8.8 -9.9 36.2   9.7 -11.4 32.9  8.9 -14.8 31.5  11.6 -17.1 30.0 9.5 -20.2 29.5  5.8 -21.0 29.3       4.3 -24.5 29.1  0.6 -25.3 29.3  0.5 -28.5 31.4  -3.2 -28.9 32.0 -6.5 -27.9 30.5 -9.5 -28.7 32.7 -9.1 -32.2 34.0 -6.7 -33.3 31.2      -3.0 -32.9 32.1 -0.4 -33.0 29.4 -0.1 -36.8 28.8 2.9 -38.3 30.7  4.2 -40.6 28.0  3.5 -38.2 25.2  6.4 -36.4 23.5  6.5 -32.9 22.0       4.6 -32.8 18.6  2.5 -35.9 19.2  -1.1 -36.5 20.3 -3.4 -39.1 21.9 -7.1 -39.7 22.0 -9.2 -39.8 25.2 -11.5 -42.8 25.8    -15.2 -42.7 26.7 -14.8 -42.3 30.4        -12.3 -39.5 30.0        -14.2 -37.7 27.3        -17.3 -38.0 29.6        -15.4 -36.5 32.6    -14.3 -33.6 30.4 -17.9 -33.1 29.4        -19.2 -33.2 33.0        -16.7 -30.4 33.8        -18.7 -28.1 31.6        -22.2 -29.4 32.6    -22.6 -31.6 29.5 -24.2 -35.0 30.3        -22.3 -38.2 29.6        -22.5 -39.9 26.3        -22.8 -43.7 25.4

1XXIJ
334
M       R       W       Y       P       W       L       R       P       D       F       E       K       L       V       A       S   L
.       .       .       .       T       T       H       H       H       H       H       H       H       H       H       T       T   .
65.7 -16.7 128.5        64.5 -19.9 130.3        66.9 -22.7 131.5        66.0 -25.8 129.7        67.1 -29.4 130.9        69.2 -30.4 127.8     71.4 -27.4 128.6        72.9 -29.3 131.5        75.3 -32.0 130.1        76.9 -29.2 128.1        76.9 -26.7 131.0        78.6 -29.4 133.1     81.8 -29.7 131.1        81.8 -26.0 130.4        82.2 -25.2 134.1        84.3 -28.3 134.6        87.0 -26.5 132.5    86.5 -23.1 134.2 87.2 -24.8 137.5        90.3 -26.9 136.4        91.9 -24.1 134.3        91.2 -25.8 130.9        88.8 -22.9 130.1    88.9 -21.4 126.6 89.0 -17.6 126.7        86.9 -17.5 123.5        83.8 -19.6 123.0        81.6 -19.9 119.9        78.4 -21.8 120.5    76.4 -23.0 117.5 72.8 -23.9 117.9        69.4 -23.9 116.3        66.9 -21.3 117.2        65.1 -21.8 120.5        67.7 -24.4 121.7    68.4 -22.5 124.8 71.4 -20.9 123.2        70.2 -17.9 125.1        70.3 -19.4 128.5        73.6 -21.3 128.3        75.7 -18.2 127.8    74.1 -16.3 130.7 74.1 -19.4 133.0        77.8 -19.8 132.3        78.8 -16.2 133.0        77.0 -16.2 136.4        79.4 -18.9 137.5    82.4 -16.9 136.6 81.1 -13.6 138.0        80.6 -14.9 141.5        82.7 -13.6 144.4        81.9 -16.6 146.7        80.7 -19.6 144.5    79.2 -22.6 146.2 80.7 -25.3 143.9        78.7 -25.6 140.7        76.6 -22.5 141.0        76.6 -19.1 142.8        75.6 -18.0 146.3    73.1 -15.2 145.4 74.5 -12.7 148.0        78.0 -11.6 146.9        78.5 -8.3 145.0 78.7 -10.2 141.8        75.6 -12.3 142.3        73.7 -9.2 143.5      74.8 -6.8 140.8 74.2 -9.8 138.5 70.6 -10.2 139.7        70.4 -6.4 139.2 71.7 -7.1 135.7 74.0 -4.0 135.8 77.6 -5.1 136.4      81.1 -3.7 135.7 82.4 -7.0 134.2 79.5 -8.2 132.2 79.4 -6.7 128.7 76.5 -7.3 126.2 77.5 -6.8 122.6 74.1 -6.4 121.0 74.1 -5.1 117.4      73.4 -1.3 116.8 69.7 -0.7 117.7 67.9 -0.2 114.4 70.3 -1.8 112.0 70.7 -5.6 111.3 74.4 -5.6 111.8 77.1 -4.7 114.4 79.7 -1.9 114.3      83.5 -1.7 114.3 82.9 1.0 117.0  81.1 -1.0 119.7 83.9 -3.6 119.5 86.8 -1.2 119.9 85.1 0.5 122.9  84.7 -2.8 124.7 88.3 -3.9 123.7      89.8 -0.9 125.5 87.2 -1.3 128.2 88.8 -4.6 129.1 92.3 -3.2 128.8 91.9 -1.8 132.3 91.6 -3.7 135.7 88.3 -5.5 136.6 86.3 -3.2 138.9      85.3 -6.1 141.2 88.7 -7.8 141.0 88.2 -11.1 139.2        87.1 -12.3 135.8        85.6 -10.1 132.9        83.1 -11.6 130.4    81.8 -10.7 126.9 78.9 -12.0 124.9        77.8 -11.6 121.3        74.2 -12.0 120.1        75.2 -13.4 116.7        78.7 -13.1 115.2    77.5 -13.8 111.7 75.7 -10.4 111.7        78.8 -8.2 112.1 80.6 -5.9 109.6 83.7 -7.5 107.9 85.8 -4.8 109.5 84.2 -5.5 112.9 85.6 -9.1 112.6      89.2 -8.1 112.3 88.8 -5.4 114.9 87.8 -8.1 117.3 90.6 -10.6 116.4        93.6 -8.6 117.7 92.2 -8.7 121.2 91.9 -12.5 121.2    95.5 -12.5 120.1 96.5 -10.2 123.0        94.2 -11.8 125.7        93.9 -9.5 128.8 94.7 -11.9 131.7        91.7 -13.9 133.1        89.1 -12.4 130.8     86.3 -14.5 129.1        84.5 -14.4 125.7        81.3 -15.9 124.4        80.5 -15.9 120.7        77.3 -17.4 119.5    76.0 -18.2 116.0 73.5 -20.3 114.2        75.8 -20.6 111.2        79.5 -20.8 112.0        80.6 -21.2 108.4        79.5 -17.8 107.3    82.1 -16.6 109.9 85.7 -15.6 109.1        87.9 -18.7 109.1        90.7 -16.8 110.8        88.2 -15.8 113.5        86.5 -19.2 114.1    90.0 -20.6 114.6 90.7 -18.5 117.7        87.6 -19.7 119.5        86.5 -23.1 120.6        83.4 -24.2 118.7        80.9 -26.1 120.8    77.9 -27.9 119.3 74.6 -27.5 121.1        73.0 -30.7 119.6        69.2 -30.8 118.8        68.0 -33.4 121.2        66.8 -36.7 119.7    63.3 -36.2 118.3 61.6 -38.6 120.7        62.8 -36.6 123.6        61.9 -33.2 122.0        58.4 -34.3 121.2        57.6 -35.6 124.7    59.0 -32.4 126.2 56.8 -30.5 123.8        53.6 -32.5 124.5        53.8 -30.9 127.9        54.2 -27.3 126.9        51.2 -27.5 124.5    48.6 -30.1 123.9 48.5 -30.9 120.2        48.7 -34.1 118.1        52.2 -35.6 118.2        52.1 -35.5 114.4        51.8 -31.8 114.5    54.3 -31.6 117.3 56.9 -33.4 115.1        56.1 -31.1 112.2        56.9 -28.0 114.3        60.2 -29.7 115.5        61.2 -30.5 111.9    60.2 -27.0 110.8 62.8 -25.7 113.3        65.7 -27.9 112.4        65.1 -29.8 115.6        65.1 -26.7 117.9        62.5 -27.5 120.5    62.4 -23.9 121.6 61.6 -22.5 118.2        59.2 -25.4 117.7        57.4 -24.4 120.9        57.2 -20.8 119.7        55.7 -22.0 116.5    52.5 -22.6 118.4 52.3 -19.7 120.9        51.0 -17.2 118.3        48.2 -18.8 116.5        50.9 -19.8 114.0        50.0 -23.4 113.8    46.4 -22.5 113.1 47.8 -20.1 110.4        49.7 -23.1 109.1        46.6 -25.3 109.0        45.0 -22.6 106.8        47.8 -23.0 104.4    47.6 -26.8 104.0 43.8 -26.5 103.4        44.3 -24.0 100.7        47.1 -26.1 99.2 44.9 -29.2 99.2 42.3 -27.4 97.2 44.7 -26.0 94.7 46.9 -29.1 94.1      44.1 -31.1 92.5 42.2 -28.1 91.1 45.1 -26.2 89.4 43.8 -22.9 90.9 47.0 -22.0 92.7 46.0 -18.4 92.8 43.6 -19.4 95.6 46.6 -19.4 98.0      46.7 -15.7 97.7 44.6 -15.2 100.8        46.6 -17.8 102.7        50.0 -16.3 101.8        49.4 -12.6 101.4        50.1 -11.5 104.9     53.0 -9.8 106.8 54.0 -13.4 107.8        54.5 -14.0 104.1        58.1 -15.2 104.4        57.5 -17.6 107.2        54.6 -18.9 105.2     56.4 -19.2 101.9        59.0 -21.3 103.7        56.2 -23.6 105.0        55.0 -24.0 101.4        58.6 -24.8 100.4    59.0 -27.6 102.9 56.0 -29.3 101.4        57.5 -28.8 97.9 60.7 -30.3 99.0 59.0 -33.4 100.5        57.0 -34.1 97.3 60.3 -33.9 95.5 62.0 -36.3 97.9      59.2 -38.8 97.5 59.9 -39.0 93.7 63.3 -40.3 94.6 61.9 -42.6 97.2 63.3 -40.7 100.2        62.2 -42.2 103.5        62.3 -39.0 105.6     59.6 -36.3 105.3        59.3 -33.1 107.2        55.6 -32.7 106.8        54.4 -35.9 108.5        51.5 -33.8 109.6    49.8 -34.3 106.2 51.3 -37.2 104.2        48.2 -37.2 102.1        49.2 -33.8 100.8        52.8 -34.5 99.9 51.5 -37.6 98.1 48.9 -36.0 95.9      51.3 -33.1 95.3 54.0 -35.5 94.1 51.3 -37.3 92.0 50.2 -34.2 90.0 53.6 -32.6 89.7 56.6 -34.0 87.9 59.9 -33.8 89.7 60.9 -31.0 87.3      57.8 -28.9 87.8 58.3 -29.2 91.5 62.0 -28.0 91.2 60.8 -25.0 89.2 58.2 -24.0 91.7 60.5 -24.7 94.7 62.8 -22.0 93.4 60.4 -19.2 92.5      58.9 -19.8 95.9 62.4 -19.6 97.2 63.7 -16.5 95.4 60.4 -14.7 95.9 60.3 -15.3 99.6 63.7 -13.7 99.7 62.8 -10.5 97.9 60.1 -10.1 100.5     62.5 -10.4 103.3        65.4 -8.5 101.8 64.0 -5.6 99.5  62.8 -2.7 101.7 59.2 -2.6 100.6 57.9 -4.3 97.4  54.2 -5.0 96.9  53.5 -8.2 98.9       50.4 -8.5 96.9  51.8 -8.1 93.4  54.8 -10.2 93.8 52.9 -13.1 95.2 50.1 -12.8 92.6 52.9 -13.0 90.1 54.5 -16.0 91.9 51.1 -17.7 91.9      50.0 -17.1 88.3 53.5 -18.0 87.3 53.3 -21.3 89.1 50.3 -22.1 87.0 52.1 -21.6 83.7 54.8 -24.0 84.7 52.0 -26.5 85.2 50.8 -26.1 81.6      52.0 -28.7 79.0 54.8 -27.5 76.7 56.2 -24.8 78.9 60.0 -24.6 79.4 60.6 -24.0 83.1 63.6 -22.0 84.5 67.2 -23.0 85.4 69.4 -21.4 88.1      71.8 -18.6 87.1 75.6 -17.8 87.3

```

Example command:

```
awk 'NR%6 != 0 {print} NR%6==0 {res=(sqrt(($1-x)^2 + ($2-y)^2 + ($3-z)^2)); if(res <= 4.5){printf "%f\n", res} else{print}} {x=$1; y=$2; z=$3}'

4LGTD
247
M       S       E       K       L       Q       K       V       L       A       R       A       G       H       G       S       R   T
.       .       E       E       H       H       H       H       H       H       H       T       T       S       S       .       H   .
-34.3 -28.9 16.3        -33.0 -26.8 13.4        -29.4 -25.4 13.3        -27.9 -22.6 11.2        -26.9 -24.1 7.8 -23.3 -22.9 8.2 -23.0 -24.7 11.5     -24.1 -27.9 9.9 -21.8 -27.7 6.9 -18.8 -26.7 9.0 -19.5 -29.7 11.3        -19.6 -31.9 8.2 -16.2 -30.7 7.2 -16.7 -28.2 4.4      -15.0 -25.1 5.8 -15.4 -22.7 8.7 -18.8 -21.3 9.7 -17.7 -17.9 8.3 -16.3 -19.3 5.1 -19.5 -21.2 4.5 -21.5 -18.0 5.1 -19.7 -16.2 2.3      -20.8 -19.1 -0.0        -24.4 -19.2 1.2 -24.6 -15.4 0.6 -23.5 -16.0 -3.0        -26.4 -18.4 -3.6        -24.1 -21.4 -4.0    -25.9 -23.6 -1.5 -29.3 -25.2 -2.0        -31.4 -27.0 0.5 -34.0 -29.5 -0.8        -33.9 -27.8 -4.2        -34.1 -24.2 -3.0        -31.3 -21.7 -2.9     -30.7 -21.0 0.8 -30.2 -17.5 2.2 -27.9 -16.3 4.9 -29.4 -17.1 8.3 -31.3 -20.1 7.0 -31.7 -23.1 9.4 -32.1 -26.8 8.7 -32.6 -30.1 10.5     -30.2 -32.9 9.7 -32.0 -36.0 8.4 -31.6 -38.7 5.7 -34.0 -36.6 3.6 -32.0 -33.4 3.7 -30.4 -32.7 0.4 -27.7 -30.0 0.6 -25.8 -29.0 -2.5     -23.0 -26.6 -3.0        -22.5 -25.6 -6.6        -25.0 -28.3 -7.5        -22.9 -31.1 -5.9        -24.6 -33.1 -3.2    -23.0 -33.1 0.2  -23.1 -36.0 2.6 -23.7 -34.9 6.2 -23.8 -36.9 9.4 -27.2 -37.3 11.1        -28.5 -37.1 14.7        -19.2 -38.7 17.8    -17.7 -36.6 20.6 -14.0 -37.3 20.3        -11.2 -35.4 22.0        -7.4 -35.3 21.2 -4.7 -34.3 23.6 -1.9 -32.6 21.6 1.6 -31.6 22.7  2.2 -28.4 20.7       5.8 -27.5 20.1  7.1 -23.9 19.9  7.1 -22.6 16.3  3.6 -24.0 15.4  0.7 -21.6 14.7  -2.8 -22.2 16.0 -4.8 -21.4 12.8 -6.4 -23.3 9.9       -4.8 -20.9 7.5  -1.6 -18.8 8.1  -1.1 -16.0 5.5  2.7 -16.4 6.0   3.0 -20.1 5.2   4.2 -21.3 8.7   3.5 -24.7 10.3  0.1 -24.9 11.9       -1.1 -27.3 14.7 -3.6 -28.9 12.3 -1.0 -30.0 9.8  -0.7 -33.4 11.5 -4.0 -34.7 12.9 -6.1 -37.9 12.3 -8.5 -37.7 9.3  -12.2 -38.1 10.1     -15.3 -39.0 8.2 -17.5 -36.0 7.6 -16.2 -34.0 10.6        -13.2 -31.8 11.6        -10.9 -31.0 14.5        -11.7 -27.8 16.4     -8.8 -26.3 18.5 -9.8 -25.6 22.1 -8.4 -22.1 22.1 -5.0 -20.9 20.9 -1.5 -21.1 22.2 1.2 -18.6 21.2  3.8 -19.7 18.6  6.6 -20.0 21.1       4.5 -21.3 23.9  4.7 -25.1 24.5  2.3 -27.7 25.8  -1.3 -29.0 25.9 -3.8 -28.1 23.2 -7.0 -30.0 23.1 -8.9 -30.7 19.9 -12.5 -32.0 19.6     -14.9 -32.8 16.9        -17.6 -31.5 19.2        -17.9 -27.7 18.8        -20.2 -27.3 21.8        -17.6 -29.0 24.0        -15.0 -26.8 22.4     -16.8 -23.5 23.1        -17.7 -24.6 26.6        -14.1 -25.3 27.4        -12.9 -21.9 26.0        -15.6 -19.6 27.2    -14.8 -18.2 30.6 -18.4 -18.5 31.9        -17.5 -22.1 32.5        -14.3 -21.2 34.4        -12.2 -24.1 33.1        -8.7 -24.1 34.5 -5.7 -23.6 32.3      -2.5 -24.5 34.2 0.9 -23.2 33.2  4.4 -23.9 34.3  6.1 -20.6 34.0  9.8 -20.0 34.3  10.1 -16.4 35.4 13.0 -14.0 35.5 13.3 -11.2 38.1      12.7 -10.5 41.8 9.7 -12.2 43.3  8.4 -12.0 46.8  5.2 -13.0 48.6  3.7 -9.5 48.4   3.8 -9.8 44.6   2.2 -13.2 44.8  -0.6 -11.7 47.0      -1.0 -8.9 44.5  -1.3 -11.3 41.5 -3.9 -13.4 43.4 -6.1 -10.5 44.4 -6.1 -8.6 41.1  -3.9 -6.2 39.3  -4.5 -3.7 36.4  -3.6 -4.8 32.9       -3.8 -2.6 29.8  -7.0 -4.4 29.0  -8.5 -4.4 32.5  -8.2 -6.2 35.8  -6.8 -9.7 35.9  -5.7 -12.2 38.6 -4.2 -15.6 38.7 -6.6 -18.1 40.2      -3.8 -20.0 41.8 -0.0 -19.4 42.0 2.4 -21.9 43.4  6.2 -21.7 43.7  7.6 -24.7 41.8  11.4 -24.5 42.2 14.2 -22.0 42.2 17.4 -20.5 41.0      18.6 -17.0 41.7 16.7 -14.6 39.4 15.1 -17.6 37.8 12.2 -19.4 39.4 9.4 -21.8 38.3  5.8 -21.1 39.2  2.7 -23.1 38.4  -0.3 -20.8 37.8      -3.8 -21.1 36.7 -6.3 -18.9 34.9 -9.9 -19.4 34.1 -9.9 -17.1 31.1 -7.6 -16.7 28.1 -6.8 -13.1 27.2 -4.0 -12.5 24.7 -0.6 -13.1 26.3      -2.2 -12.9 29.7 0.3 -14.9 31.8  3.2 -12.8 30.5  1.3 -9.5 31.1   0.4 -10.5 34.6  4.0 -11.2 35.6  5.1 -8.0 33.8   2.5 -6.0 35.8        4.6 -6.8 38.9   8.0 -6.3 37.2   8.8 -9.9 36.2   9.7 -11.4 32.9  8.9 -14.8 31.5  11.6 -17.1 30.0 9.5 -20.2 29.5  5.8 -21.0 29.3       4.3 -24.5 29.1  0.6 -25.3 29.3  0.5 -28.5 31.4  -3.2 -28.9 32.0 -6.5 -27.9 30.5 -9.5 -28.7 32.7 -9.1 -32.2 34.0 -6.7 -33.3 31.2      -3.0 -32.9 32.1 -0.4 -33.0 29.4 -0.1 -36.8 28.8 2.9 -38.3 30.7  4.2 -40.6 28.0  3.5 -38.2 25.2  6.4 -36.4 23.5  6.5 -32.9 22.0       4.6 -32.8 18.6  2.5 -35.9 19.2  -1.1 -36.5 20.3 -3.4 -39.1 21.9 -7.1 -39.7 22.0 -9.2 -39.8 25.2 -11.5 -42.8 25.8    -15.2 -42.7 26.7 -14.8 -42.3 30.4        -12.3 -39.5 30.0        -14.2 -37.7 27.3        -17.3 -38.0 29.6        -15.4 -36.5 32.6    -14.3 -33.6 30.4 -17.9 -33.1 29.4        -19.2 -33.2 33.0        -16.7 -30.4 33.8        -18.7 -28.1 31.6        -22.2 -29.4 32.6    -22.6 -31.6 29.5 -24.2 -35.0 30.3        -22.3 -38.2 29.6        -22.5 -39.9 26.3        -22.8 -43.7 25.4

1XXIJ
334
M       R       W       Y       P       W       L       R       P       D       F       E       K       L       V       A       S   L
.       .       .       .       T       T       H       H       H       H       H       H       H       H       H       T       T   .
65.7 -16.7 128.5        64.5 -19.9 130.3        66.9 -22.7 131.5        66.0 -25.8 129.7        67.1 -29.4 130.9        69.2 -30.4 127.8     71.4 -27.4 128.6        72.9 -29.3 131.5        75.3 -32.0 130.1        76.9 -29.2 128.1        76.9 -26.7 131.0        78.6 -29.4 133.1     81.8 -29.7 131.1        81.8 -26.0 130.4        82.2 -25.2 134.1        84.3 -28.3 134.6        87.0 -26.5 132.5    86.5 -23.1 134.2 87.2 -24.8 137.5        90.3 -26.9 136.4        91.9 -24.1 134.3        91.2 -25.8 130.9        88.8 -22.9 130.1    88.9 -21.4 126.6 89.0 -17.6 126.7        86.9 -17.5 123.5        83.8 -19.6 123.0        81.6 -19.9 119.9        78.4 -21.8 120.5    76.4 -23.0 117.5 72.8 -23.9 117.9        69.4 -23.9 116.3        66.9 -21.3 117.2        65.1 -21.8 120.5        67.7 -24.4 121.7    68.4 -22.5 124.8 71.4 -20.9 123.2        70.2 -17.9 125.1        70.3 -19.4 128.5        73.6 -21.3 128.3        75.7 -18.2 127.8    74.1 -16.3 130.7 74.1 -19.4 133.0        77.8 -19.8 132.3        78.8 -16.2 133.0        77.0 -16.2 136.4        79.4 -18.9 137.5    82.4 -16.9 136.6 81.1 -13.6 138.0        80.6 -14.9 141.5        82.7 -13.6 144.4        81.9 -16.6 146.7        80.7 -19.6 144.5    79.2 -22.6 146.2 80.7 -25.3 143.9        78.7 -25.6 140.7        76.6 -22.5 141.0        76.6 -19.1 142.8        75.6 -18.0 146.3    73.1 -15.2 145.4 74.5 -12.7 148.0        78.0 -11.6 146.9        78.5 -8.3 145.0 78.7 -10.2 141.8        75.6 -12.3 142.3        73.7 -9.2 143.5      74.8 -6.8 140.8 74.2 -9.8 138.5 70.6 -10.2 139.7        70.4 -6.4 139.2 71.7 -7.1 135.7 74.0 -4.0 135.8 77.6 -5.1 136.4      81.1 -3.7 135.7 82.4 -7.0 134.2 79.5 -8.2 132.2 79.4 -6.7 128.7 76.5 -7.3 126.2 77.5 -6.8 122.6 74.1 -6.4 121.0 74.1 -5.1 117.4      73.4 -1.3 116.8 69.7 -0.7 117.7 67.9 -0.2 114.4 70.3 -1.8 112.0 70.7 -5.6 111.3 74.4 -5.6 111.8 77.1 -4.7 114.4 79.7 -1.9 114.3      83.5 -1.7 114.3 82.9 1.0 117.0  81.1 -1.0 119.7 83.9 -3.6 119.5 86.8 -1.2 119.9 85.1 0.5 122.9  84.7 -2.8 124.7 88.3 -3.9 123.7      89.8 -0.9 125.5 87.2 -1.3 128.2 88.8 -4.6 129.1 92.3 -3.2 128.8 91.9 -1.8 132.3 91.6 -3.7 135.7 88.3 -5.5 136.6 86.3 -3.2 138.9      85.3 -6.1 141.2 88.7 -7.8 141.0 88.2 -11.1 139.2        87.1 -12.3 135.8        85.6 -10.1 132.9        83.1 -11.6 130.4    81.8 -10.7 126.9 78.9 -12.0 124.9        77.8 -11.6 121.3        74.2 -12.0 120.1        75.2 -13.4 116.7        78.7 -13.1 115.2    77.5 -13.8 111.7 75.7 -10.4 111.7        78.8 -8.2 112.1 80.6 -5.9 109.6 83.7 -7.5 107.9 85.8 -4.8 109.5 84.2 -5.5 112.9 85.6 -9.1 112.6      89.2 -8.1 112.3 88.8 -5.4 114.9 87.8 -8.1 117.3 90.6 -10.6 116.4        93.6 -8.6 117.7 92.2 -8.7 121.2 91.9 -12.5 121.2    95.5 -12.5 120.1 96.5 -10.2 123.0        94.2 -11.8 125.7        93.9 -9.5 128.8 94.7 -11.9 131.7        91.7 -13.9 133.1        89.1 -12.4 130.8     86.3 -14.5 129.1        84.5 -14.4 125.7        81.3 -15.9 124.4        80.5 -15.9 120.7        77.3 -17.4 119.5    76.0 -18.2 116.0 73.5 -20.3 114.2        75.8 -20.6 111.2        79.5 -20.8 112.0        80.6 -21.2 108.4        79.5 -17.8 107.3    82.1 -16.6 109.9 85.7 -15.6 109.1        87.9 -18.7 109.1        90.7 -16.8 110.8        88.2 -15.8 113.5        86.5 -19.2 114.1    90.0 -20.6 114.6 90.7 -18.5 117.7        87.6 -19.7 119.5        86.5 -23.1 120.6        83.4 -24.2 118.7        80.9 -26.1 120.8    77.9 -27.9 119.3 74.6 -27.5 121.1        73.0 -30.7 119.6        69.2 -30.8 118.8        68.0 -33.4 121.2        66.8 -36.7 119.7    63.3 -36.2 118.3 61.6 -38.6 120.7        62.8 -36.6 123.6        61.9 -33.2 122.0        58.4 -34.3 121.2        57.6 -35.6 124.7    59.0 -32.4 126.2 56.8 -30.5 123.8        53.6 -32.5 124.5        53.8 -30.9 127.9        54.2 -27.3 126.9        51.2 -27.5 124.5    48.6 -30.1 123.9 48.5 -30.9 120.2        48.7 -34.1 118.1        52.2 -35.6 118.2        52.1 -35.5 114.4        51.8 -31.8 114.5    54.3 -31.6 117.3 56.9 -33.4 115.1        56.1 -31.1 112.2        56.9 -28.0 114.3        60.2 -29.7 115.5        61.2 -30.5 111.9    60.2 -27.0 110.8 62.8 -25.7 113.3        65.7 -27.9 112.4        65.1 -29.8 115.6        65.1 -26.7 117.9        62.5 -27.5 120.5    62.4 -23.9 121.6 61.6 -22.5 118.2        59.2 -25.4 117.7        57.4 -24.4 120.9        57.2 -20.8 119.7        55.7 -22.0 116.5    52.5 -22.6 118.4 52.3 -19.7 120.9        51.0 -17.2 118.3        48.2 -18.8 116.5        50.9 -19.8 114.0        50.0 -23.4 113.8    46.4 -22.5 113.1 47.8 -20.1 110.4        49.7 -23.1 109.1        46.6 -25.3 109.0        45.0 -22.6 106.8        47.8 -23.0 104.4    47.6 -26.8 104.0 43.8 -26.5 103.4        44.3 -24.0 100.7        47.1 -26.1 99.2 44.9 -29.2 99.2 42.3 -27.4 97.2 44.7 -26.0 94.7 46.9 -29.1 94.1      44.1 -31.1 92.5 42.2 -28.1 91.1 45.1 -26.2 89.4 43.8 -22.9 90.9 47.0 -22.0 92.7 46.0 -18.4 92.8 43.6 -19.4 95.6 46.6 -19.4 98.0      46.7 -15.7 97.7 44.6 -15.2 100.8        46.6 -17.8 102.7        50.0 -16.3 101.8        49.4 -12.6 101.4        50.1 -11.5 104.9     53.0 -9.8 106.8 54.0 -13.4 107.8        54.5 -14.0 104.1        58.1 -15.2 104.4        57.5 -17.6 107.2        54.6 -18.9 105.2     56.4 -19.2 101.9        59.0 -21.3 103.7        56.2 -23.6 105.0        55.0 -24.0 101.4        58.6 -24.8 100.4    59.0 -27.6 102.9 56.0 -29.3 101.4        57.5 -28.8 97.9 60.7 -30.3 99.0 59.0 -33.4 100.5        57.0 -34.1 97.3 60.3 -33.9 95.5 62.0 -36.3 97.9      59.2 -38.8 97.5 59.9 -39.0 93.7 63.3 -40.3 94.6 61.9 -42.6 97.2 63.3 -40.7 100.2        62.2 -42.2 103.5        62.3 -39.0 105.6     59.6 -36.3 105.3        59.3 -33.1 107.2        55.6 -32.7 106.8        54.4 -35.9 108.5        51.5 -33.8 109.6    49.8 -34.3 106.2 51.3 -37.2 104.2        48.2 -37.2 102.1        49.2 -33.8 100.8        52.8 -34.5 99.9 51.5 -37.6 98.1 48.9 -36.0 95.9      51.3 -33.1 95.3 54.0 -35.5 94.1 51.3 -37.3 92.0 50.2 -34.2 90.0 53.6 -32.6 89.7 56.6 -34.0 87.9 59.9 -33.8 89.7 60.9 -31.0 87.3      57.8 -28.9 87.8 58.3 -29.2 91.5 62.0 -28.0 91.2 60.8 -25.0 89.2 58.2 -24.0 91.7 60.5 -24.7 94.7 62.8 -22.0 93.4 60.4 -19.2 92.5      58.9 -19.8 95.9 62.4 -19.6 97.2 63.7 -16.5 95.4 60.4 -14.7 95.9 60.3 -15.3 99.6 63.7 -13.7 99.7 62.8 -10.5 97.9 60.1 -10.1 100.5     62.5 -10.4 103.3        65.4 -8.5 101.8 64.0 -5.6 99.5  62.8 -2.7 101.7 59.2 -2.6 100.6 57.9 -4.3 97.4  54.2 -5.0 96.9  53.5 -8.2 98.9       50.4 -8.5 96.9  51.8 -8.1 93.4  54.8 -10.2 93.8 52.9 -13.1 95.2 50.1 -12.8 92.6 52.9 -13.0 90.1 54.5 -16.0 91.9 51.1 -17.7 91.9      50.0 -17.1 88.3 53.5 -18.0 87.3 53.3 -21.3 89.1 50.3 -22.1 87.0 52.1 -21.6 83.7 54.8 -24.0 84.7 52.0 -26.5 85.2 50.8 -26.1 81.6      52.0 -28.7 79.0 54.8 -27.5 76.7 56.2 -24.8 78.9 60.0 -24.6 79.4 60.6 -24.0 83.1 63.6 -22.0 84.5 67.2 -23.0 85.4 69.4 -21.4 88.1      71.8 -18.6 87.1 75.6 -17.8 87.3

```

NB. in both of these 'blocks' the equation is \>4.5, so the entire line is printed. Is this your desired outcome?
> 2 votes

---
Tags: protein-structure, linux, awk
---


# Question
Title: GLM for multinomial using r - results of genotyping
I have next data:

* age - how old is person, e.g. 17 or 25
* bmi - body mass index, e.g. 18.25 or 29.35
* rsXXX - results of genotyping in rs(point), e.g. rs10852521: C/C, C/T and T/T

I devide all data in 2 groups `bmi_group`, create new column bmi\_group, remove all rows not in group

* low - bmi\<25.0
* high - bmi\>=30.0

```
LG = 25.0
HG = 30.0
subdata = subset(all_data, bmi < LG | bmi >= HG)
subdata$bmi_group <- with(subdata, ifelse(bmi < LG, 'low', 'high'))

```

I need to find out if a certain allelic variant affects the risk of obesity, taking into account age. I try next:

```
> model <- glm(imt_group ~ rs11075990+age, data = point.data, family = binomial)
> summary(model)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-2.4209  -0.6968   0.4239   0.6470   2.0662

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)    6.07127    0.56333  10.778   <2e-16 ***
rs11075990A/G -0.65093    0.28982  -2.246   0.0247 *
rs11075990G/G -0.53818    0.35860  -1.501   0.1334
age           -0.12783    0.01235 -10.351   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 568.58  on 430  degrees of freedom
Residual deviance: 403.53  on 427  degrees of freedom
AIC: 411.53

Number of Fisher Scoring iterations: 4

```

Problem is genotyping data not binominal - `Factor w/ 3 levels "A/A","A/G","G/G":`. What and how can I do?

# Answer
I do not see a problem here. There are 3 levels corresponding to the three genotypes and genotype `A/A` is the base (reference) level and the "effects" of the other two genotypes are stated "on top of the base level".

This tutorial on linear regression with factor levels is very useful to understand how to conduct linear regression in general but also goes over the specifics of doing so with R.

Moreover, you created a new discrete (factor) variable out of a continuous variable when you bin BMI status as low or high. You can certainly do so but you can also use it as it is, as a continuous variable. And I guess you forgot to include this variable in your model.
> 3 votes

---
Tags: r, statistics, molecular-genetics
---