# Question
Title: Multiple genome assemblies of the same bacterial species
I have some RNA-seq data where there are reads from the "host" as well as from several bacteria species. In this experimental context, I am interested in the host associated reads and the latter is regarded as "contaminants".

I am using `BBSplit` to split host- and bacteria-associated reads using my genome of interest as well as a "repo" containing thousands of bacterial genomes that are likely to be present in the host.

This repo has multiple genome assemblies for the same species making the `BBSplit` step computationally expensive. The total size of the repo is approximatly 25GB, more than 25X the human genome! I am wondering what would be the best criteria to pick one assembly of a particular genome over another assembly.

All in all, different genome assemblies of a given species have more or less similar total lengths, %GC, number of contigs meaning that picking any would probably be fine but not being a microbiologist, I would like to hear expert's opinions on this. My current take is to take the assembly with the largest size for each species.

# Answer
A couple of ideas:

1. Use kraken2 to classify the reads. If you are able to convert your "repo" to kraken2 format, then the matching will be very fast and relatively accurate. Also, using bwa-mem on a joint database could work.
2. I am not a big fan of custom and not well defined "repos". Why don't you just use the official RefSeq repository of prokaryote genomes? (https://www.ncbi.nlm.nih.gov/refseq/about/prokaryotes/). I don't remember its size, but I remember using it without any problem.
> 1 votes

# Answer
1. What I would suggest, without knowing much biology about the system is it's a tree problem. The difficulty is knowing whether the repo is prokaryotic. To solve it I'd take the 16S gene and make a tree - each 16S sequence contains its genome ID. I'd then draw a tree, e.g. IQ-TREE is good and fast. This will describe the overall diversity and provides and informed basis of how to represent each species.
2. Your approach certainly has merit because again if these are bacteria they go through gene loss - which can cause problems. Selecting the largest genome ensures you minimise this loss.

What I would do is combine 1 and 2. This is to label the taxa that has the largest genome and then produce the 16S tree. It will then be obvious by manual inspection if your approach is representative of the genetic diversity. If it isn't .... you can do the whole thing manually (could be painful), the other approach is to use a tree reading package/module/library. I use ETE3, you'd use ape (@haci's from memory is a big R codes). The issue is traversing the tree from the 'biggest genome' to the most distant genetically related member of the species. What you do is identify the MRCA (most recent common ancestor) for that species and pick a taxa that is in the opposite clade to the biggest genome and any member will do. Thus this is automating the process via Python or R code so thousands of genomes/species can be screened.

To try and be clearer, what you are doing in the code is moving the tree one node at a time and checking if all the members below it are the same species. As soon as the answer is 'no', it's the previous node thats the MRCA for the species. In a bifurcating tree that divides the data into its two biggest groups of genetic diversity.

This approach ensures you get a reasonable representation of the species, by selecting two members of it. One is the largest genome.

It's a while since I've used ETE3, so I can't remember the methods, there's probably a MRCA method for a given species. You can do this empirically, doing a recursive climb up the tree from the biggest genome until your monophyly contains a difference species. The previous traversal is the MRCA, then you ask it to divide the species into two clades and any member of the clade without the biggest clade is what you want.

I agree that it needs a bit of "tree-thinking" to do this ... but it's doable. You then use the genome IDs to filter your repo and produce a greatly slimmed down repo.

You could alternative switch to an HPC method (e.g. AWS) using the entire repo

The alternative approach is know a lot about the biology/ pathogenesis of the problem.
> -1 votes

---
Tags: rna-seq, assembly, microbiology
---


# Question
Title: How do I quantify a specific somatic variant?
I'm working with targeted Illumina sequencing data generated with DNA from diseased and healthy tissue (this is an age-related disease and is not cancer/neoplastic). My hypothesis is that the diseased tissue might contain specific somatic mutations in one particular gene (let's call it GENE), while the healthy tissue will not. The reasoning behind this: The same mutations cause autosomal-dominant versions of my disease of interest when present in the germline (my diseased tissue is from sporadic disease). My coverage varies, but is in the range of 1000x - 50'000x. My sequencing protocol employed UMIs.

I have run Mutect2 in default settings and with lower thresholds (`--tumor-lod-to-emit`) and haven't found any calls that correspond to these pathogenic mutations, so my results are negative. Nevertheless, I have reason to believe that my pathogenic mutations might be present at **very** low VAF (around 0.05-0.1 %), based on other experiments.

**My question**: How do I quantify the number of reads that correspond to a particular variant of interest? Let's say that base 2000 in GENE is usually an A, whereas the mutation A\>T would be pathogenic. How do I count the number of Ts in position 2000?

As I'm somewhat new to bioinformatics, I'm also very open to any other ideas. Thanks!

---

**Edit**: @terdon's approach worked perfectly in cases where I used a VCF for `--alleles` that had been generated as an output of Mutect2, such as the following one:

*VCF produced by Mutect*:

```
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT
chr20   4699605 .   A   G   .   .   AS_SB_TABLE=881,932|920,847;DP=3791;ECNT=1;MBQ=20,17;MFRL=248,243;MMQ=60,60;MPOS=21;POPAF=7.30;TLOD=2014.57 GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:1813,1767:0.469:3580:559,218:384,191:1720,1679:881,932,920,847

```

Curiously, if I used a manually produced file such as the following, no output was produced, unless the variant of interest was heterozygous:

*Manually produced VCF*:

```
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT
chr20   4699605 rs1799990   A   G   .   .   .

```

*Output using the manual VCF, heterozygous sample*:

```
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  heterogyzous_sample
chr20   4699605 .   A   G   .   .   AS_SB_TABLE=881,932|920,847;DP=3791;ECNT=1;MBQ=20,17;MFRL=248,243;MMQ=60,60;MPOS=21;POPAF=7.30;TLOD=2014.57 GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:1813,1767:0.469:3580:559,218:384,191:1720,1679:881,932,920,847

```

*Output using the manual VCF, non-heterozygous sample*:

```
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  novariant_sample

```

Manual inspection of this non-heterozygous sample revealed that it contained 1 read with the ALT variant, so this doesn't seem to depend on the variant being absent.

# Answer
I have a little function that can report this by using samtools to look at the read pileup. Um. It isn't the most elegant piece of code, but it does sorta do the job:

```
variantSupport ()
{
  genome=$1
  bam=$2;
  range=$3;
  samtools mpileup -f "$genome" -r "$range" "$bam" | cut -f 5 |
    tr '[a-z]' '[A-Z]' | fold -w 1 | sort | uniq -c |
    awk 'BEGIN{tot=1}$2~/,/{ $2="." }
         {
           tot+=$1;
           c[$2]+=$1
         }
         $2~/[a-zA-Z]/{ varReads+=$1; }
         END{
           for(i in c){
             print i,c[i]
           }
           printf "%s of %s reads (%.1f%%), support the variant.\n",varReads,tot,varReads*100/tot
   }'
}

```

Add the function to your `~/.bashrc` or just paste it directly into a terminal and you can then do:

```
$ variantSupport /path/to/genome/hg19.fa foo.bam chr13:28607916-28607916
[mpileup] 1 samples in 1 input files
A 915
* 1
. 93
915 of 1010 reads (90.6%), support the variant.

```

This is telling me that at position 28607916 of chr13, I have 1 read supporting a deletion (`*`, see https://www.htslib.org/doc/samtools-mpileup.html for an explanation of the symbols), 93 reads supporting the reference residue (T, in this case) and 915 supporting an A, so a T=\> A variant.

However, a better approach is to use a variant caller and force it to call the variant you want. You can then investigate the support for it. With `mutect2`, you can use the `--alleles` option to force the calling and then use the `-L` option to limit calling to only the specific region. If you give the same vcf file for both options, that will return results for the variants you want and only those. For example, given this file `foo.vcf.gz` (which needs to be compressed using `bgzip` and indexed using `tabix`):

```
$ zcat 1.vcf.gz
##fileformat=VCFv4.3
##reference=/genomes/hg19/hg19.fa
##FORMAT=<ID=AB,Number=A,Type=Float,Description="Allelic Balance">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read depth">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Net Genotype quality across all datasets, calculated from GQ scores of callsets supporting the consensus GT, using only one callset from each dataset">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr13,assembly=hg19,length=115169878>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  HG003
chr13   28607917    .   G   C   .   .   .   GT:AB:DP:GQ 1/1:0.6:395:318

```

You can then do:

```
gatk Mutect2 \
  -R /genomes/hg19/hg19.fa \
  --alleles 1.vcf.gz  \
  -L 1.vcf.gz \
  -I foo.bam \
  -O out.vcf

```

This will produce an `out.vcf` like this (I am excluding most headers for clarity):

```

##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  331970
chr13   28607917    .   G   C   .   .   DP=1107;ECNT=1;MBQ=30,0;MFRL=353,0;MMQ=60,60;MPOS=50;POPAF=7.30;TLOD=-2.716e+00 GT:AD:AF:DP:F1R2:F2R1:SB    0/1:1043,0:9.616e-04:1043:518,0:496,0:989,54,0,0

```

I can see this is not a real variant since the value for AD, the number of reads supporting the reference and the variant is `1043,0`, meaning that 1043 reads support the reference and none support the variant.

This way, you can get results for your variants of interest and see if there is any evidence supporting them.
> 4 votes

---
Tags: variant-calling, illumina, gatk
---


# Question
Title: Calculating Fisher's exact test for COG categories
This is a continuation of the following question -

Fisher's exact test for COG categories in pan, core genome analysis

Also, the dataset -

How do I perform a Fisher's exact test on this data to show that a COG category is significantly more or less in the pan/core/accessory/species specific genomes?

# Answer
What I would personally do is run a Χ2

```
from scipy.stats import chisquare
import numpy as np

pan = [11, 39, 4 ....]
core = [29, 50, 7 ... ]
f_exp = np.array(pan)
f_obs = np. array(core)
print (chisquare(f_obs=f_obs, f_exp=f_exp))

```

This will output the total statistic and whether its above or below the critical probability.

Fisher's Exact test is essentially Χ2 for a small number of observations, which is not the case with your data.

Looking at your data Χ2 is not likely to be "significant" (which is the result you are looking for), i.e. the pan genome is the core genome across all COGs.
> 1 votes

---
Tags: statistics, eggnog, pangenome
---


# Question
Title: Why is mtDNA not labelled as such for many organisms?
Also posted on biostars

While I was searching the genomes of some eukariotic organisms (fungi,plants,protists) even though I knew that they had mitochondrial DNA,when I searched their genome structure I couldn't find any "chromosome" labeled as "mtDNA". Here I give some links of genomes that show this: https://www.ncbi.nlm.nih.gov/genome/gdv/browser/genome/?id=GCF\_016861865.1 https://www.ncbi.nlm.nih.gov/genome/gdv/browser/genome/?id=GCF\_009556855.1

Why does this happens? How can I find if the mtDNA sequence is inluded in these "chromosomes"?

# Answer
The most obvious answer is wet-lab rather than bioinformatics: mtDNA can get removed in the DNA extraction process. As a control group, in bacterial drug-resistance studies the plasmids are frequently lost. mtDNA at ~15KB is smaller than a regular bacterial plasmid.

On the bioinformatics front there would be an issue about reference genome assemblies as well, i.e. if mtDNA genome isn't present then it will not appear as a separate genome.
> 1 votes

---
Tags: genome, sequence-annotation, reference-genome, mtdna
---


title: All against all sequence alignment

text:
# Question
Title: All against all sequence alignment
I have a large set of protein sequences and would like to create a "similarity matrix" out of them (for phylo trees later). So I have to run `n(n-1)` sequence alignments. Wondering what library I should use for that and also what sophisticated lib would provide this all against all without me writing a for loop and parsing results myself...

Currently using `Bio.pairwise2.align.globalxx` from BioPython as a development placeholder but I do not think the results look really useful? It returns a `score > 1` and no similarity in `[0,1[` or `evalue`.

I have been out of the Bioinfo game for a while, so any suggestions for libraries (python preferred) much appreciated! :)

# Answer
> 1 votes
I wouldn't use Python here. I've looked at using `Bio.pairwise2.align.` and its has major shortcomings. In this context the problem is establishing homologous alignment positions ... if you want to make a tree.

Whats the problems? This is you start introducing the gap character '-' as an additional state in constructing your tree and the tree will bias towards not the maximum number of substitutions but the maximum gap (indel) difference between the sequences.

The second problem is changes in homologous sites between different pairwise alignments. This gets complicated fast.

Is that bad? The problem is we have little understanding of the mutation model behind generating indels and there isn't a model to account for varying descriptions of homologous alignment positions. However, we have a very clear understanding of amino acid substitutions in a single global alignment, i.e. a singular description of homologous amino acid positions. A lot of work has gone into modelling amino acid mutation, with a tree a practical output of that work.

---

Thus, I would simply align using muscle5 under the **-super5** option. This will produce a global alignment very quickly and for *vast* amounts of data. This alignment will produce the best model of homologous amino acid sites.

If you just want a distance matrix then import the entire alignment into MEGAXI and request the JTT model of amino acid substitution with the gamma distribution. This will output a distance matrix nicely and will then allow you to produce a tree.

BTW estimating the gamma distribution without maximum likelihood can be difficult in distance methods (MEGA might have resolved this). This tree will be much better than the one you are proposing (you can estimate this parameter via e.g. IQTREE or ModelTest).

Why? Again alot of work has been done on modelling amino acid substations and JTT is one of the better models (assuming this is nuclear sequence not mtDNA). The LG matrix is much better BTW, but I don't think thats in MEGA.

---
Tags: sequence-alignment, multiple-sequence-alignment, biopython
---


# Question
Title: How can I download multiple SRR* files at a time when using sratoolkit's fastq-dump?
I have multiple SRR files to download, and their accession is in a list:

```
SRR9004325
SRR9004326
SRR9004327

```

I read the NCBI document but don't find how to achieve this using the command "fastq-dump".

And I noticed that if I use the command `prefetch`, I can download from the list file:

```
prefetch --option-file SraAccList.txt

```

So is there any parameter in "fastq-dump" to download accessions stored in a file?

As @terdon said, it is better to use `fasterq-dump`, and the command: `fasterq-dump --option-file SraAccList.txt` doesn't work. When I try it, I just get the usage message:

```
$ fasterq-dump --option-file SraAccList.txt

Usage:
  fasterq-dump <path> [options]
  fasterq-dump <accession> [options]

Options:
  -F|--format                      format (special, fastq, default=fastq)
  -o|--outfile                     output-file
  -O|--outdir                      output-dir
  -b|--bufsize                     size of file-buffer dflt=1MB
  -c|--curcache                    size of cursor-cache dflt=10MB
  -m|--mem                         memory limit for sorting dflt=100MB
  -t|--temp                        where to put temp. files dflt=curr dir
  -e|--threads                     how many thread dflt=6
  -p|--progress                    show progress
  -x|--details                     print details
  -s|--split-spot                  split spots into reads
  -S|--split-files                 write reads into different files
  -3|--split-3                     writes single reads in special file
  --concatenate-reads              writes whole spots into one file
  -Z|--stdout                      print output to stdout
  -f|--force                       force to overwrite existing file(s)
  --skip-technical                 skip technical reads
  --include-technical              include technical reads
  -M|--min-read-len                filter by sequence-len
  --table                          which seq-table to use in case of pacbio
  -B|--bases                       filter by bases
  -A|--append                      append to output-file
  --fasta                          produce FASTA output
  --fasta-unsorted                 produce FASTA output, unsorted
  --fasta-ref-tbl                  produce FASTA output from REFERENCE tbl
  --fasta-concat-all               concatenate all rows and produce FASTA
  --internal-ref                   extract only internal REFERENCEs
  --external-ref                   extract only external REFERENCEs
  --ref-name                       extract only these REFERENCEs
  --ref-report                     enumerate references
  --use-name                       print name instead of seq-id
  --seq-defline                    custom defline for sequence:  $ac=accession,
                               $sn=spot-name,  $sg=spot-group, $si=spot-id,
                                   $ri=read-id, $rl=read-length
  --qual-defline                   custom defline for qualities:  same as
                                   seq-defline
  -U|--only-unaligned              process only unaligned reads
  -a|--only-aligned                process only aligned reads
  --disk-limit                     explicitly set disk-limit
  --disk-limit-tmp                 explicitly set disk-limit for temp. files
  --size-check                     switch to control: on=perform size-check
                                   (default),  off=do not perform size-check,
                                   only=perform size-check only
  --ngc <PATH>                     PATH to ngc file

  -h|--help                        Output brief explanation for the program.
  -V|--version                     Display the version of the program then
                                   quit.
  -L|--log-level <level>           Logging level as number or enum string. One
                                   of (fatal|sys|int|err|warn|info|debug) or
                                   (0-6) Current/default is warn.
  -v|--verbose                     Increase the verbosity of the program
                                   status messages. Use multiple times for more
                                   verbosity. Negates quiet.
  -q|--quiet                       Turn off all status messages for the
                                   program. Negated by verbose.
  --option-file <file>             Read more options and parameters from the
                                   file.
for more information visit:
   https://github.com/ncbi/sra-tools/wiki/HowTo:-fasterq-dump
   https://github.com/ncbi/sra-tools/wiki/08.-prefetch-and-fasterq-dump

fasterq-dump : 3.0.6

```

# Answer
> 2 votes
If `fasterq-dump` can't read from a file, you can simply loop over the file yourself:

```
while read srr; do
  fasterq-dump "$srr" -o "$srr".fastq
done < SraAccList.txt

```

That will iterate over your `SraAccList.txt` file and extract the fastqs for each accession. I tested with this command:

```
fasterq-dump SRR9004325 -o ff.fastq

```

And that created files `ff_1.fastq` and `ff_2.fastq`, so you can expect files named like `SRR9004325_1.fastq` and `SRR9004325_1.fastq` if you use the loop above.

You could try and parallelize it but I doubt it is worth it since the bottleneck will always be the download speed.

---
Tags: data-download, ncbi
---


title: How do I summarise VCF files?

text:
# Question
Title: How do I summarise VCF files?
I have got a number of VCF files that all contain data on one particular variant, with every VCF being created from an individual sample. So there is one entry per file (headers omitted):

```
##source=Mutect2
##tumor_sample=312091_5-Ctrl1_G01_S6_R1_001.fastq.gz
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  312091_5-Ctrl1_G01_S6_R1_001.fastq.gz
chr20   4699818 .   G   A   .   .   AS_SB_TABLE=805,864|0,0;DP=1871;ECNT=1;MBQ=17,0;MFRL=273,0;MMQ=60,60;MPOS=50;POPAF=7.30;TLOD=-3.537e+00 GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:1669,0:6.794e-04:1669:252,0:233,0:1595,0:805,864,0,0

```

I would like to create a summary table of the parameters given within the VCFs, including of the parameters within the `FORMAT` field. Ideally this table would have the following as header:

```
Sample CHROM    POS ID  REF ALT QUAL    FILTER  GT INFO-DP FORMAT-DP AD-REF AD-ALT

```

and every sample/variant combination would be a line. I would like to use such a table to create statistics for my samples (on zygosity, number of ALT variants etc).

Is there a pre-made tool that is capable of doing this? I'm also open to any other suggestions. Thanks!

**Edit**: Here are two VCF files (I have about 30 in total) and the desired output

**First VCF**:

```
##fileformat=VCFv4.2
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=AF,Number=A,Type=Float,Description="Allele fractions of alternate alleles in the tumor">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##FORMAT=<ID=F1R2,Number=R,Type=Integer,Description="Count of reads in F1R2 pair orientation supporting each allele">
##FORMAT=<ID=F2R1,Number=R,Type=Integer,Description="Count of reads in F2R1 pair orientation supporting each allele">
##FORMAT=<ID=FAD,Number=R,Type=Integer,Description="Count of fragments supporting each allele.">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=PGT,Number=1,Type=String,Description="Physical phasing haplotype information, describing how the alternate alleles are phased in relation to one another; will always be heterozygous and is not intended to describe called alleles">
##FORMAT=<ID=PID,Number=1,Type=String,Description="Physical phasing ID information, where each unique ID within a given sample (but not across samples) connects records within a phasing group">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
##FORMAT=<ID=PS,Number=1,Type=Integer,Description="Phasing set (typically the position of the first variant in the set)">
##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias.">
##GATKCommandLine=<ID=Mutect2,CommandLine="Mutect2 --alleles /add/2024-03-04_force_mutect/E200K.vcf.gz --output /add/2024-03-04_force_mutect/output//Ctrl5_E200K_frequency.vcf --intervals /add/2024-03-04_force_mutect/E200K.vcf.gz --input /add/2023-06-18_CJD_8_samples/312091_8-Ctrl5_H02_S3_R1_001.fastq.gz.picard.markedDup.recal.bam --reference /home/mcarta/databases/hg38.fa --f1r2-median-mq 50 --f1r2-min-bq 20 --f1r2-max-depth 200 --flow-likelihood-parallel-threads 0 --flow-likelihood-optimized-comp false --trim-to-haplotype true --exact-matching false --flow-use-t0-tag false --flow-probability-threshold 0.003 --flow-remove-non-single-base-pair-indels false --flow-remove-one-zero-probs false --flow-quantization-bins 121 --flow-fill-empty-bins-value 0.001 --flow-symmetric-indel-probs false --flow-report-insertion-or-deletion false --flow-disallow-probs-larger-than-call false --flow-lump-probs false --flow-retain-max-n-probs-base-format false --flow-probability-scaling-factor 10 --flow-order-cycle-length 4 --keep-boundary-flows false --genotype-pon-sites false --genotype-germline-sites false --af-of-alleles-not-in-resource -1.0 --mitochondria-mode false --mutect3-training-mode false --mutect3-ref-downsample 10 --mutect3-alt-downsample 20 --mutect3-non-artifact-ratio 20 --tumor-lod-to-emit 3.0 --initial-tumor-lod 2.0 --pcr-snv-qual 40 --pcr-indel-qual 40 --base-qual-correction-factor 5 --max-population-af 0.01 --downsampling-stride 1 --callable-depth 10 --max-suspicious-reads-per-alignment-start 0 --normal-lod 2.2 --ignore-itr-artifacts false --gvcf-lod-band -2.5 --gvcf-lod-band -2.0 --gvcf-lod-band -1.5 --gvcf-lod-band -1.0 --gvcf-lod-band -0.5 --gvcf-lod-band 0.0 --gvcf-lod-band 0.5 --gvcf-lod-band 1.0 --minimum-allele-fraction 0.0 --independent-mates false --flow-mode NONE --disable-adaptive-pruning false --kmer-size 10 --kmer-size 25 --dont-increase-kmer-sizes-for-cycles false --allow-non-unique-kmers-in-ref false --num-pruning-samples 1 --min-dangling-branch-length 4 --recover-all-dangling-branches false --max-num-haplotypes-in-population 128 --min-pruning 2 --adaptive-pruning-initial-error-rate 0.001 --pruning-lod-threshold 2.302585092994046 --pruning-seeding-lod-threshold 9.210340371976184 --max-unpruned-variants 100 --linked-de-bruijn-graph false --disable-artificial-haplotype-recovery false --enable-legacy-graph-cycle-detection false --debug-assembly false --debug-graph-transformations false --capture-assembly-failure-bam false --num-matching-bases-in-dangling-end-to-recover -1 --error-correction-log-odds -Infinity --error-correct-reads false --kmer-length-for-read-error-correction 25 --min-observations-for-kmer-to-be-solid 20 --likelihood-calculation-engine PairHMM --base-quality-score-threshold 18 --dragstr-het-hom-ratio 2 --dont-use-dragstr-pair-hmm-scores false --pair-hmm-gap-continuation-penalty 10 --expected-mismatch-rate-for-read-disqualification 0.02 --pair-hmm-implementation FASTEST_AVAILABLE --pcr-indel-model CONSERVATIVE --phred-scaled-global-read-mismapping-rate 45 --disable-symmetric-hmm-normalizing false --disable-cap-base-qualities-to-map-quality false --enable-dynamic-read-disqualification-for-genotyping false --dynamic-read-disqualification-threshold 1.0 --native-pair-hmm-threads 4 --native-pair-hmm-use-double-precision false --flow-hmm-engine-min-indel-adjust 6 --flow-hmm-engine-flat-insertion-penatly 45 --flow-hmm-engine-flat-deletion-penatly 45 --pileup-detection false --use-pdhmm false --use-pdhmm-overlap-optimization false --make-determined-haps-from-pd-code false --print-pileupcalling-status false --fallback-gga-if-pdhmm-fails true --pileup-detection-enable-indel-pileup-calling false --pileup-detection-active-region-phred-threshold 0.0 --num-artificial-haplotypes-to-add-per-allele 5 --artifical-haplotype-filtering-kmer-size 10 --pileup-detection-snp-alt-threshold 0.1 --pileup-detection-indel-alt-threshold 0.1 --pileup-detection-absolute-alt-depth 0.0 --pileup-detection-snp-adjacent-to-assembled-indel-range 5 --pileup-detection-snp-basequality-filter 12 --pileup-detection-bad-read-tolerance 0.0 --pileup-detection-proper-pair-read-badness true --pileup-detection-edit-distance-read-badness-threshold 0.08 --pileup-detection-chimeric-read-badness true --pileup-detection-template-mean-badness-threshold 0.0 --pileup-detection-template-std-badness-threshold 0.0 --pileup-detection-filter-assembly-alt-bad-read-tolerance 0.0 --pileup-detection-edit-distance-read-badness-for-assembly-filtering-threshold 0.12 --bam-writer-type CALLED_HAPLOTYPES --dont-use-soft-clipped-bases false --override-fragment-softclip-check false --min-base-quality-score 10 --smith-waterman FASTEST_AVAILABLE --emit-ref-confidence NONE --max-mnp-distance 1 --force-call-filtered-alleles false --reference-model-deletion-quality 30 --soft-clip-low-quality-ends false --allele-informative-reads-overlap-margin 2 --smith-waterman-dangling-end-match-value 25 --smith-waterman-dangling-end-mismatch-penalty -50 --smith-waterman-dangling-end-gap-open-penalty -110 --smith-waterman-dangling-end-gap-extend-penalty -6 --smith-waterman-haplotype-to-reference-match-value 200 --smith-waterman-haplotype-to-reference-mismatch-penalty -150 --smith-waterman-haplotype-to-reference-gap-open-penalty -260 --smith-waterman-haplotype-to-reference-gap-extend-penalty -11 --smith-waterman-read-to-haplotype-match-value 10 --smith-waterman-read-to-haplotype-mismatch-penalty -15 --smith-waterman-read-to-haplotype-gap-open-penalty -30 --smith-waterman-read-to-haplotype-gap-extend-penalty -5 --flow-assembly-collapse-hmer-size 0 --flow-assembly-collapse-partial-mode false --flow-filter-alleles false --flow-filter-alleles-qual-threshold 30.0 --flow-filter-alleles-sor-threshold 3.0 --flow-filter-lone-alleles false --flow-filter-alleles-debug-graphs false --min-assembly-region-size 50 --max-assembly-region-size 300 --active-probability-threshold 0.002 --max-prob-propagation-distance 50 --force-active false --assembly-region-padding 100 --padding-around-indels 75 --padding-around-snps 20 --padding-around-strs 75 --max-extension-into-assembly-region-padding-legacy 25 --max-reads-per-alignment-start 50 --enable-legacy-assembly-region-trimming false --interval-set-rule UNION --interval-padding 0 --interval-exclusion-padding 0 --interval-merging-rule ALL --read-validation-stringency SILENT --seconds-between-progress-updates 10.0 --disable-sequence-dictionary-validation false --create-output-bam-index true --create-output-bam-md5 false --create-output-variant-index true --create-output-variant-md5 false --max-variants-per-shard 0 --lenient false --add-output-sam-program-record true --add-output-vcf-command-line true --cloud-prefetch-buffer 40 --cloud-index-prefetch-buffer -1 --disable-bam-index-caching false --sites-only-vcf-output false --help false --version false --showHidden false --verbosity INFO --QUIET false --use-jdk-deflater false --use-jdk-inflater false --gcs-max-retries 20 --gcs-project-for-requester-pays  --disable-tool-default-read-filters false --max-read-length 2147483647 --min-read-length 30 --minimum-mapping-quality 20 --disable-tool-default-annotations false --enable-all-annotations false",Version="4.5.0.0",Date="March 5, 2024 at 10:20:03 AM UTC">
##INFO=<ID=AS_SB_TABLE,Number=1,Type=String,Description="Allele-specific forward/reverse read counts for strand bias tests. Includes the reference and alleles separated by |.">
##INFO=<ID=AS_UNIQ_ALT_READ_COUNT,Number=A,Type=Integer,Description="Number of reads with unique start and mate end positions for each alt at a variant site">
##INFO=<ID=CONTQ,Number=1,Type=Float,Description="Phred-scaled qualities that alt allele are not due to contamination">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth; some reads may have been filtered">
##INFO=<ID=ECNT,Number=1,Type=Integer,Description="Number of events in this haplotype">
##INFO=<ID=GERMQ,Number=1,Type=Integer,Description="Phred-scaled quality that alt alleles are not germline variants">
##INFO=<ID=MBQ,Number=R,Type=Integer,Description="median base quality by allele">
##INFO=<ID=MFRL,Number=R,Type=Integer,Description="median fragment length by allele">
##INFO=<ID=MMQ,Number=R,Type=Integer,Description="median mapping quality by allele">
##INFO=<ID=MPOS,Number=A,Type=Integer,Description="median distance from end of read">
##INFO=<ID=NALOD,Number=A,Type=Float,Description="Negative log 10 odds of artifact in normal with same allele fraction as tumor">
##INFO=<ID=NCount,Number=1,Type=Integer,Description="Count of N bases in the pileup">
##INFO=<ID=NLOD,Number=A,Type=Float,Description="Normal log 10 likelihood ratio of diploid het or hom alt genotypes">
##INFO=<ID=OCM,Number=1,Type=Integer,Description="Number of alt reads whose original alignment doesn't match the current contig.">
##INFO=<ID=PON,Number=0,Type=Flag,Description="site found in panel of normals">
##INFO=<ID=POPAF,Number=A,Type=Float,Description="negative log 10 population allele frequencies of alt alleles">
##INFO=<ID=ROQ,Number=1,Type=Float,Description="Phred-scaled qualities that alt allele are not due to read orientation artifact">
##INFO=<ID=RPA,Number=R,Type=Integer,Description="Number of times tandem repeat unit is repeated, for each allele (including reference)">
##INFO=<ID=RU,Number=1,Type=String,Description="Tandem repeat unit (bases)">
##INFO=<ID=SEQQ,Number=1,Type=Integer,Description="Phred-scaled quality that alt alleles are not sequencing errors">
##INFO=<ID=STR,Number=0,Type=Flag,Description="Variant is a short tandem repeat">
##INFO=<ID=STRANDQ,Number=1,Type=Integer,Description="Phred-scaled quality of strand bias artifact">
##INFO=<ID=STRQ,Number=1,Type=Integer,Description="Phred-scaled quality that alt alleles in STRs are not polymerase slippage errors">
##INFO=<ID=TLOD,Number=A,Type=Float,Description="Log 10 likelihood ratio score of variant existing versus not existing">
##MutectVersion=2.2
##contig=<ID=chr2,length=242193529>
##contig=<ID=chr20,length=64444167>
##contig=<ID=chr4,length=190214555>
##filtering_status=Warning: unfiltered Mutect 2 calls.  Please run FilterMutectCalls to remove false positives.
##source=Mutect2
##tumor_sample=312091_8-Ctrl5_H02_S3_R1_001.fastq.gz
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  312091_8-Ctrl5_H02_S3_R1_001.fastq.gz
chr20   4699818 .   G   A   .   .   AS_SB_TABLE=1376,1402|0,0;DP=2939;ECNT=1;MBQ=18,0;MFRL=257,0;MMQ=60,60;MPOS=50;POPAF=7.30;TLOD=-3.455e+00   GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:2778,0:4.100e-04:2778:472,0:398,0:2493,0:1376,1402,0,0

```

**Second VCF**:

```
##fileformat=VCFv4.2
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=AF,Number=A,Type=Float,Description="Allele fractions of alternate alleles in the tumor">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##FORMAT=<ID=F1R2,Number=R,Type=Integer,Description="Count of reads in F1R2 pair orientation supporting each allele">
##FORMAT=<ID=F2R1,Number=R,Type=Integer,Description="Count of reads in F2R1 pair orientation supporting each allele">
##FORMAT=<ID=FAD,Number=R,Type=Integer,Description="Count of fragments supporting each allele.">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=PGT,Number=1,Type=String,Description="Physical phasing haplotype information, describing how the alternate alleles are phased in relation to one another; will always be heterozygous and is not intended to describe called alleles">
##FORMAT=<ID=PID,Number=1,Type=String,Description="Physical phasing ID information, where each unique ID within a given sample (but not across samples) connects records within a phasing group">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
##FORMAT=<ID=PS,Number=1,Type=Integer,Description="Phasing set (typically the position of the first variant in the set)">
##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias.">
##GATKCommandLine=<ID=Mutect2,CommandLine="Mutect2 --alleles /add/2024-03-04_force_mutect/E200K.vcf.gz --output /add/2024-03-04_force_mutect/output//Ctrl4_E200K_frequency.vcf --intervals /add/2024-03-04_force_mutect/E200K.vcf.gz --input /add/2023-06-18_CJD_16_samples/312092_15-Ctrl4_H01_S11_R1_001.fastq.gz.picard.markedDup.recal.bam --reference /home/mcarta/databases/hg38.fa --f1r2-median-mq 50 --f1r2-min-bq 20 --f1r2-max-depth 200 --flow-likelihood-parallel-threads 0 --flow-likelihood-optimized-comp false --trim-to-haplotype true --exact-matching false --flow-use-t0-tag false --flow-probability-threshold 0.003 --flow-remove-non-single-base-pair-indels false --flow-remove-one-zero-probs false --flow-quantization-bins 121 --flow-fill-empty-bins-value 0.001 --flow-symmetric-indel-probs false --flow-report-insertion-or-deletion false --flow-disallow-probs-larger-than-call false --flow-lump-probs false --flow-retain-max-n-probs-base-format false --flow-probability-scaling-factor 10 --flow-order-cycle-length 4 --keep-boundary-flows false --genotype-pon-sites false --genotype-germline-sites false --af-of-alleles-not-in-resource -1.0 --mitochondria-mode false --mutect3-training-mode false --mutect3-ref-downsample 10 --mutect3-alt-downsample 20 --mutect3-non-artifact-ratio 20 --tumor-lod-to-emit 3.0 --initial-tumor-lod 2.0 --pcr-snv-qual 40 --pcr-indel-qual 40 --base-qual-correction-factor 5 --max-population-af 0.01 --downsampling-stride 1 --callable-depth 10 --max-suspicious-reads-per-alignment-start 0 --normal-lod 2.2 --ignore-itr-artifacts false --gvcf-lod-band -2.5 --gvcf-lod-band -2.0 --gvcf-lod-band -1.5 --gvcf-lod-band -1.0 --gvcf-lod-band -0.5 --gvcf-lod-band 0.0 --gvcf-lod-band 0.5 --gvcf-lod-band 1.0 --minimum-allele-fraction 0.0 --independent-mates false --flow-mode NONE --disable-adaptive-pruning false --kmer-size 10 --kmer-size 25 --dont-increase-kmer-sizes-for-cycles false --allow-non-unique-kmers-in-ref false --num-pruning-samples 1 --min-dangling-branch-length 4 --recover-all-dangling-branches false --max-num-haplotypes-in-population 128 --min-pruning 2 --adaptive-pruning-initial-error-rate 0.001 --pruning-lod-threshold 2.302585092994046 --pruning-seeding-lod-threshold 9.210340371976184 --max-unpruned-variants 100 --linked-de-bruijn-graph false --disable-artificial-haplotype-recovery false --enable-legacy-graph-cycle-detection false --debug-assembly false --debug-graph-transformations false --capture-assembly-failure-bam false --num-matching-bases-in-dangling-end-to-recover -1 --error-correction-log-odds -Infinity --error-correct-reads false --kmer-length-for-read-error-correction 25 --min-observations-for-kmer-to-be-solid 20 --likelihood-calculation-engine PairHMM --base-quality-score-threshold 18 --dragstr-het-hom-ratio 2 --dont-use-dragstr-pair-hmm-scores false --pair-hmm-gap-continuation-penalty 10 --expected-mismatch-rate-for-read-disqualification 0.02 --pair-hmm-implementation FASTEST_AVAILABLE --pcr-indel-model CONSERVATIVE --phred-scaled-global-read-mismapping-rate 45 --disable-symmetric-hmm-normalizing false --disable-cap-base-qualities-to-map-quality false --enable-dynamic-read-disqualification-for-genotyping false --dynamic-read-disqualification-threshold 1.0 --native-pair-hmm-threads 4 --native-pair-hmm-use-double-precision false --flow-hmm-engine-min-indel-adjust 6 --flow-hmm-engine-flat-insertion-penatly 45 --flow-hmm-engine-flat-deletion-penatly 45 --pileup-detection false --use-pdhmm false --use-pdhmm-overlap-optimization false --make-determined-haps-from-pd-code false --print-pileupcalling-status false --fallback-gga-if-pdhmm-fails true --pileup-detection-enable-indel-pileup-calling false --pileup-detection-active-region-phred-threshold 0.0 --num-artificial-haplotypes-to-add-per-allele 5 --artifical-haplotype-filtering-kmer-size 10 --pileup-detection-snp-alt-threshold 0.1 --pileup-detection-indel-alt-threshold 0.1 --pileup-detection-absolute-alt-depth 0.0 --pileup-detection-snp-adjacent-to-assembled-indel-range 5 --pileup-detection-snp-basequality-filter 12 --pileup-detection-bad-read-tolerance 0.0 --pileup-detection-proper-pair-read-badness true --pileup-detection-edit-distance-read-badness-threshold 0.08 --pileup-detection-chimeric-read-badness true --pileup-detection-template-mean-badness-threshold 0.0 --pileup-detection-template-std-badness-threshold 0.0 --pileup-detection-filter-assembly-alt-bad-read-tolerance 0.0 --pileup-detection-edit-distance-read-badness-for-assembly-filtering-threshold 0.12 --bam-writer-type CALLED_HAPLOTYPES --dont-use-soft-clipped-bases false --override-fragment-softclip-check false --min-base-quality-score 10 --smith-waterman FASTEST_AVAILABLE --emit-ref-confidence NONE --max-mnp-distance 1 --force-call-filtered-alleles false --reference-model-deletion-quality 30 --soft-clip-low-quality-ends false --allele-informative-reads-overlap-margin 2 --smith-waterman-dangling-end-match-value 25 --smith-waterman-dangling-end-mismatch-penalty -50 --smith-waterman-dangling-end-gap-open-penalty -110 --smith-waterman-dangling-end-gap-extend-penalty -6 --smith-waterman-haplotype-to-reference-match-value 200 --smith-waterman-haplotype-to-reference-mismatch-penalty -150 --smith-waterman-haplotype-to-reference-gap-open-penalty -260 --smith-waterman-haplotype-to-reference-gap-extend-penalty -11 --smith-waterman-read-to-haplotype-match-value 10 --smith-waterman-read-to-haplotype-mismatch-penalty -15 --smith-waterman-read-to-haplotype-gap-open-penalty -30 --smith-waterman-read-to-haplotype-gap-extend-penalty -5 --flow-assembly-collapse-hmer-size 0 --flow-assembly-collapse-partial-mode false --flow-filter-alleles false --flow-filter-alleles-qual-threshold 30.0 --flow-filter-alleles-sor-threshold 3.0 --flow-filter-lone-alleles false --flow-filter-alleles-debug-graphs false --min-assembly-region-size 50 --max-assembly-region-size 300 --active-probability-threshold 0.002 --max-prob-propagation-distance 50 --force-active false --assembly-region-padding 100 --padding-around-indels 75 --padding-around-snps 20 --padding-around-strs 75 --max-extension-into-assembly-region-padding-legacy 25 --max-reads-per-alignment-start 50 --enable-legacy-assembly-region-trimming false --interval-set-rule UNION --interval-padding 0 --interval-exclusion-padding 0 --interval-merging-rule ALL --read-validation-stringency SILENT --seconds-between-progress-updates 10.0 --disable-sequence-dictionary-validation false --create-output-bam-index true --create-output-bam-md5 false --create-output-variant-index true --create-output-variant-md5 false --max-variants-per-shard 0 --lenient false --add-output-sam-program-record true --add-output-vcf-command-line true --cloud-prefetch-buffer 40 --cloud-index-prefetch-buffer -1 --disable-bam-index-caching false --sites-only-vcf-output false --help false --version false --showHidden false --verbosity INFO --QUIET false --use-jdk-deflater false --use-jdk-inflater false --gcs-max-retries 20 --gcs-project-for-requester-pays  --disable-tool-default-read-filters false --max-read-length 2147483647 --min-read-length 30 --minimum-mapping-quality 20 --disable-tool-default-annotations false --enable-all-annotations false",Version="4.5.0.0",Date="March 5, 2024 at 10:21:02 AM UTC">
##INFO=<ID=AS_SB_TABLE,Number=1,Type=String,Description="Allele-specific forward/reverse read counts for strand bias tests. Includes the reference and alleles separated by |.">
##INFO=<ID=AS_UNIQ_ALT_READ_COUNT,Number=A,Type=Integer,Description="Number of reads with unique start and mate end positions for each alt at a variant site">
##INFO=<ID=CONTQ,Number=1,Type=Float,Description="Phred-scaled qualities that alt allele are not due to contamination">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth; some reads may have been filtered">
##INFO=<ID=ECNT,Number=1,Type=Integer,Description="Number of events in this haplotype">
##INFO=<ID=GERMQ,Number=1,Type=Integer,Description="Phred-scaled quality that alt alleles are not germline variants">
##INFO=<ID=MBQ,Number=R,Type=Integer,Description="median base quality by allele">
##INFO=<ID=MFRL,Number=R,Type=Integer,Description="median fragment length by allele">
##INFO=<ID=MMQ,Number=R,Type=Integer,Description="median mapping quality by allele">
##INFO=<ID=MPOS,Number=A,Type=Integer,Description="median distance from end of read">
##INFO=<ID=NALOD,Number=A,Type=Float,Description="Negative log 10 odds of artifact in normal with same allele fraction as tumor">
##INFO=<ID=NCount,Number=1,Type=Integer,Description="Count of N bases in the pileup">
##INFO=<ID=NLOD,Number=A,Type=Float,Description="Normal log 10 likelihood ratio of diploid het or hom alt genotypes">
##INFO=<ID=OCM,Number=1,Type=Integer,Description="Number of alt reads whose original alignment doesn't match the current contig.">
##INFO=<ID=PON,Number=0,Type=Flag,Description="site found in panel of normals">
##INFO=<ID=POPAF,Number=A,Type=Float,Description="negative log 10 population allele frequencies of alt alleles">
##INFO=<ID=ROQ,Number=1,Type=Float,Description="Phred-scaled qualities that alt allele are not due to read orientation artifact">
##INFO=<ID=RPA,Number=R,Type=Integer,Description="Number of times tandem repeat unit is repeated, for each allele (including reference)">
##INFO=<ID=RU,Number=1,Type=String,Description="Tandem repeat unit (bases)">
##INFO=<ID=SEQQ,Number=1,Type=Integer,Description="Phred-scaled quality that alt alleles are not sequencing errors">
##INFO=<ID=STR,Number=0,Type=Flag,Description="Variant is a short tandem repeat">
##INFO=<ID=STRANDQ,Number=1,Type=Integer,Description="Phred-scaled quality of strand bias artifact">
##INFO=<ID=STRQ,Number=1,Type=Integer,Description="Phred-scaled quality that alt alleles in STRs are not polymerase slippage errors">
##INFO=<ID=TLOD,Number=A,Type=Float,Description="Log 10 likelihood ratio score of variant existing versus not existing">
##MutectVersion=2.2
##contig=<ID=chr2,length=242193529>
##contig=<ID=chr20,length=64444167>
##contig=<ID=chr4,length=190214555>
##filtering_status=Warning: unfiltered Mutect 2 calls.  Please run FilterMutectCalls to remove false positives.
##source=Mutect2
##tumor_sample=312092_15-Ctrl4_H01_S11_R1_001.fastq.gz
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  312092_15-Ctrl4_H01_S11_R1_001.fastq.gz
chr20   4699818 .   G   A   .   .   AS_SB_TABLE=22,12|0,0;DP=48;ECNT=1;MBQ=11,0;MFRL=243,0;MMQ=60,60;MPOS=50;POPAF=7.30;TLOD=-1.584e+00 GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:34,0:0.026:34:0,0:0,0:33,0:22,12,0,0

```

**Desired output** (values from the first VCF above, with tabs as delimiter)

```
Sample CHROM    POS ID  REF ALT QUAL    FILTER  GT INFO-DP FORMAT-DP AD-REF AD-ALT
312091_8-Ctrl5_H02_S3_R1_001.fastq.gz chr20 4699818 . G A . . 0/1 2939 2778 2778 0

```

I suppose merging them would be an option. Is there a method to extract the information I'd like from a merged VCF?

# Answer
The first step is to merge your VCFs. For this, you need to compress them with `bgzip` and then index with `tabix`:

```
for file in *vcf; do
  bgzip "$file"
  tabix "$file".gz
done

```

Now, feed the files to `bcftools merge`. In this example, I have saved your two files as `sample1.vcf.gz` and `sample2.vcf.gz`, and they're the only vcf files in the directory, so I can catch them all with `*vcf.gz`:

```
bcftools merge *.vcf.gz > merged.vcf

```

This results in a file containing a single variant line, with an entry for each of the two samples:

```
chr20   4699818 .   G   A   .   .   AS_SB_TABLE=1376,1402|0,0;ECNT=1;DP=2987;MBQ=11,0;MFRL=243,0;MMQ=60,60;MPOS=50;POPAF=7.3;TLOD=-1.584    GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:2778,0:0.00041:2778:472,0:398,0:2493,0:1376,1402,0,0    0/1:34,0:0.026:34:0,0:0,0:33,0:22,12,0,0

```

From this, `bcftools query` can give you *almost* exactly what you need. The main problem is that when merging, the INFO fields are not combined. Only one is kept. In my test, it was the one from the last file processed, but this could be random, I don't know. This is the closest I could get to your desired format:

```
$ { printf 'Sample\tCHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tGT\tFORMAT-DP\tAD\n';
  bcftools query --format "[%SAMPLE ]\t%CHROM\t%POS0\t%ID\t%REF\t%ALT\t%QUAL\t%FILTER\t[%GT ]\t[%DP ]\t[%AD ] " merged.vcf ; } | column -s$'\t' -t
Sample                                                                          CHROM  POS      ID  REF  ALT  QUAL  FILTER  GT          FORMAT-DP  AD
312091_8-Ctrl5_H02_S3_R1_001.fastq.gz 312092_15-Ctrl4_H01_S11_R1_001.fastq.gz   chr20  4699817  .   G    A    .     .        0/1  0/1   2778 34    2778,0 34,0

```

The main command is this:

```
bcftools query --format "[%SAMPLE ]\t%CHROM\t%POS0\t%ID\t%REF\t%ALT\t%QUAL\t%FILTER\t[%GT ]\t[%DP ]\t[%AD ] " merged.vcf

```

The rest just add a header and pretty formatting for clarity. Please see https://samtools.github.io/bcftools/bcftools.html#query for details on the options used.

This is close, but not quite. You have the sample names listed on the same line, we're missing the DP value from the FORMAT field and we have the ADs shown in pairs, not in different columns as requested. A little post-processing can fix everything except the missing DP/FORMAT:

```
$ {
  printf 'Sample\tCHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tGT\tFORMAT-DP\tAD-REF\tAD-ALT\n';
  bcftools query --format "[%SAMPLE ]\t%CHROM\t%POS0\t%ID\t%REF\t%ALT\t%QUAL\t%FILTER\t[%GT ]\t[%DP ]\t[%AD ] " merged.vcf |
   awk -F'\t' -v OFS="\t" '
  {
    split($1, samples, " ");
    split($9, genotypes, " ");
    split($10, dps, " ");
    split($11, ads, " ");
    for(s in samples){
       split(ads[s], sampleAd, ",");
       print samples[s], $2, $3, $4, $5, $6, $7, $8, genotypes[s], dps[s], sampleAd[1], sampleAd[2]
     }
  }';
} | column -s$'\t' -t
Sample                                   CHROM  POS      ID       REF  ALT  QUAL  FILTER  GT   FORMAT-DP  AD-REF  AD-ALT
312091_8-Ctrl5_H02_S3_R1_001.fastq.gz    chr20  4699817  .  G    A    .     .       0/1  2778       2778    0
312092_15-Ctrl4_H01_S11_R1_001.fastq.gz  chr20  4699817 .  G    A    .     .       0/1  34         34      0

```

**Important**: the command above assumes there are no spaces in your sample names nor in any of the values returned. It also assumes only one ALT allele per variant. Please update your question if these assumptions are not valid.

This is as close as I can get it using "standard" tools and relatively little post-processing. It is everything you asked for except the INFO DP. If you absolutely must have that too (which seems a bit odd, honestly, you usually want the one actually reported for the variant), I fear that will require just parsing the whole thing manually.
> 2 votes

# Answer
The simplest alternative is to have a directory with all the vcf files you want to merge, and use bcftools merge. Example:

```
vcf_dir=path_to_your_vcf_directory
bcftools merge $vcf_dir/*.vcf -o bcftools_merged.vcf

```
> -1 votes

---
Tags: vcf, genomics
---


title: Creating custom VCF file

text:
# Question
Title: Creating custom VCF file
I'm stuck on what seems to be a banal problem, but I can't work out where the issue is. I need to create some customised VCF files that specify particular variants of interest.

This is a custom VCF I created for this variant:

```
##fileformat=VCFv4.3
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=20240308
##source=Ensembl
##reference=GRCh38
#CHROM  POS ID  REF ALT QUAL  FILTER    INFO
20  4699752 rs74315403  G   A   .   .  .

```

The formatting seems correct to me; the values are separated by tabs, all relevant column headers are present, there is a value for every header. I try to use it in the Mutect2 `--alleles` and `-L` options, as follows:

```
gatk Mutect2 \
  -R ref.fa \
  --alleles my.vcf.gz  \
  -L my.vcf.gz\
  -I input.bam \
  -O output.vcf

```

This gives me the following error:

```
Caused by: htsjdk.tribble.TribbleException$MalformedFeatureFile: Unable to parse header with error: Your input file has a malformed header: there are not enough columns present in the header line: #CHROM     POS     ID      REF     ALT     QUAL  FILTER    INFO, for input source: path/to/my.vcf
        at htsjdk.tribble.TribbleIndexedFeatureReader.readHeader(TribbleIndexedFeatureReader.java:265)
        at htsjdk.tribble.TribbleIndexedFeatureReader.<init>(TribbleIndexedFeatureReader.java:104)
        at htsjdk.tribble.TribbleIndexedFeatureReader.<init>(TribbleIndexedFeatureReader.java:129)
        at htsjdk.tribble.AbstractFeatureReader.getFeatureReader(AbstractFeatureReader.java:121)
        at org.broadinstitute.hellbender.engine.FeatureDataSource.getTribbleFeatureReader(FeatureDataSource.java:433)
        ... 15 more
Caused by: htsjdk.tribble.TribbleException$InvalidHeader: Your input file has a malformed header: there are not enough columns present in the header line: #CHROM       POS     ID      REF    ALT      QUAL  FILTER    INFO
        at htsjdk.variant.vcf.AbstractVCFCodec.parseHeaderFromLines(AbstractVCFCodec.java:149)
        at htsjdk.variant.vcf.VCFCodec.readActualHeader(VCFCodec.java:117)
        at htsjdk.tribble.AsciiFeatureCodec.readHeader(AsciiFeatureCodec.java:79)
        at htsjdk.tribble.AsciiFeatureCodec.readHeader(AsciiFeatureCodec.java:37)
        at htsjdk.tribble.TribbleIndexedFeatureReader.readHeader(TribbleIndexedFeatureReader.java:263)

```

Supposedly there's something wrong with the header. Can you spot the issue? Thanks!

# Answer
> 3 votes
That sounds like your header isn't tab separated. You can check for the presence of tabs with `od -c`. So, to check the header, do

```
zgrep '^#C' my.vcf.gz | od -c

```

Then check the output, it should show `\t` between each header field, like this:

```
$ zgrep '^#C'  variants.vcf.gz | od -c
0000000   #   C   H   R   O   M  \t   P   O   S  \t   I   D  \t   R   E
0000020   F  \t   A   L   T  \t   Q   U   A   L  \t   F   I   L   T   E
0000040   R  \t   I   N   F   O  \t   F   O   R   M   A   T  \t   E   D
0000060   2   0   0   1   3   9   6  \n
0000070

```

---
Tags: vcf, genomics, mutect2
---


# Question
Title: Downloading multiple SRA files as FASTQ files using SRA Toolkit
I was able to download all of the SRA files for a specific bioproject using the command:

```
prefetch -O output_directory --option-file SRR_Acc_List.txt

```

I then tried the following command to download those SRA files into FASTQ files:

```
find -name \*.sra -exec bsub fastq-dump --split-3 --gzip {} \;

```

It gave the error:

```
find: illegal option -- n
usage: find [-H | -L | -P] [-EXdsx] [-f path] path ... [expression]
       find [-H | -L | -P] [-EXdsx] -f path [path ...] [expression]

```

I then tried to add the path:

```
find /Users/lucykelly/Desktop/Assel_paper/output_directory -name '*.sra' -exec bsub fastq-dump --split-3 --gzip {} \;

```

It gave the error:

```
find: bsub: No such file or directory

```

I then tried

```
find . -name '*.sra' -exec fasterq-dump --gzip {} \;

```

and it gave the error:

```
find: fasterq-dump: No such file or directory

```

# Answer
That is a classical problem where the folder containing the binary of the tool is not in PATH, see for example https://askubuntu.com/questions/109381/how-to-add-path-of-a-program-to-path-environment-variable or many other threads towards PATH in the SE network.
> 1 votes

---
Tags: fastq, sratoolkit, sra, fastq-dump
---


title: Finding insertion in a plasmid

text:
# Question
Title: Finding insertion in a plasmid
*This question was also asked on Biostars*

I received a set of Nanopore long reads of a ~3kbp long known plasmid with an unknown insert between 1kbp and 7kbp. Could you suggest a pipeline to assemble the reads and find the insertion?

# Answer
When I've needed to do this, I use canu to assemble the nanopore reads (choosing the highest-quality reads that are at least half the length of the target), then map the assembled consensus back to the reference sequence.

The Canu assembly part of that is described here:

https://www.protocols.io/view/plasmid-sequence-analysis-from-long-reads-36wgq4n5yvk5?step=12

There are numerous different ways to map the consensus assembled sequence to the reference sequence. BLASTn in "align 2 or more sequences" mode could be used, for example. Another option, because the sequences are relatively small and you're just aligning 2, is needle from the EMBOSS suite.
> 2 votes

---
Tags: nanopore, long-reads
---


title: rna seq matrix design hybrids

text:
# Question
Title: rna seq matrix design hybrids
I am struggling to get my head around the correct design to use for my experiment. I have RNAseq of 2 pure species (s1, s2) and of the hybrid (h). I aligned all three genotypes to the s1 genomes. I am interested in the comparison of hybrids vs s1 (maternal genome). But I want to remove a potential bias from reads from s2 chromosome being harder to align to s1. (The hybrid genome is composed of both s1 and s2.) I would like to correct my DE test by removing potential effect s2 vs s1. (I aligned s2 to s1 and got the count.) But my issue is, I fail to have a full rank design. I do not have any strong feelings about the library (deseq2, edgeR...).

Following the tutorial

below is the condition matrix, as you can see if I want to make the contrast ~ s2 + condition. The matrix is not full rank. I believe I could write my design as ~s2 + h - s1 ?? But I don't know how to use deseq2 for this approach.

Any approach/tips to achieve my goal here? Is it event possible? what

| sample | condition | s2 |
| --- | --- | --- |
| 1 | s1 | 0 |
| 2 | s1 | 0 |
| 3 | s1 | 0 |
| 4 | s2 | 1 |
| 5 | s2 | 1 |
| 6 | s2 | 1 |
| 7 | h | 0 |
| 8 | h | 0 |
| 9 | h | 0 |

any help is warmly appreciated!

# Answer
I think you are over thinking the problem. You need to classify the genes of interest and to do that a classifier such as eggNOG is used. This will remove any bias. I understand eggNOG because it's based on protein phylogenetics and phylogenetics is what I do and I can assure you its very solid and the algorithm nicely constructed.

DESeq2 (I ain't a huge fan BTW, but avoid controversy) will happily accept eggNOG classification, here

Then its simply three tests

1. S1 vs S2
2. S1 vs hybrid
3. S2 vs hybrid

I personally would perform unsupervised learning, particularly via Scanpy first on S1, S2 and hybrid separately. This will give a rough indication of the key difference of expression profiles between S1 and S2 and see where the hydrid occurs.
> 1 votes

---
Tags: rna-seq, deseq2, edger, design
---


# Question
Title: Extracting a cytochrome B sequence from NCBI's nucleotide database
Can someone tell the way to extract the fasta sequences for the gene cytb of *Acetes japonicus* (shrimp important to China and South Korea)?

Can I extract them directly from NCBI nucleotide database (i.e. nuccore)?

For instance, I'm trying to fetch the fasta data of *Acetes japonicus* of cytb gene. So I'm using biopython, like this:

```
 handle = Entrez.esearch(db="nucleotide", term="Acetes japonicus[Orgn] AND cytb[Gene]")
record=Entrez.read(handle)

```

And the output doesn't give the an ID, however when I search manually I get the result.

# Answer
> 0 votes
I should have answered correctly ages back, because I do know it.

Cytochrome b is a massively used target for arthropods - these are shrimps. However, the nomenclature changed:

* Traditionally nomenclature `cytb`
* 'New' nomenclature `COB`

Wait you might ask, arn't you confusing that with COI, II, III? No because that nomenclature changed too:

* Traditional nomenclature `COI`
* 'New' nomenclature `COX1`

'New' was a long time ago though.

Thus the gene you are looking for is MZ571562, here

```
 ... AND COB [gene]

```

# Answer
> 0 votes
Here is how you can do it from command line version:

```
esearch -db nuccore -query "COB [GENE] AND "Acetes japonicus" [ORGN]"|efetch -format gb

```

OR

```
esearch -db nuccore -query "COB [GENE] AND "Acetes japonicus" [ORGN]"|efetch -format fasta

```

Note that gene name is actually COB.

---
Tags: database, nucleotide-models, efetch
---


# Question
Title: How to perform meaningful Gene set erichment analysis or otherwise find broader themes/functions in different cell populations?
(For context: I'm somewhat new to bioinformatic analyses but am mostly comfortable with R) I have identified transcriptomically and morphologically different cancer cell populations within a patient (scRNAseq, spatial transcriptomics, IHC). They are quite distinct, with hundreds of genes as differentially expressed genes (DEGs) even after high tresholds. It's quite interesting to describe, but I am looking to find better ways to represent their features. Of course I am attaching DEG tables as supplementary data etc. but I thought that gene set enrichment analysis could give broad terms that would a) help me understand which DEGs might contribute together to biological processes and b) describe the bigger picture.

So, I tried running a few terms (GO:BP, GO:CC, GO:MF, ...) from the GSEA msigdb website on the DEGs from scRNAseq, but have come out a bit frustrated. A lot of terms make sense and reflect the changes, but have large overlap. For example, the 3 top GO:BP terms (Extracellular matrix organization", "Extracellular structure organization", "External encapsulating structure organziation" are scored identically, with identical DEGs as the genes that are part of their gene sets. This is repeated further down the list with subsets of this gene set. Now, reporting these doesn't actually do what I want, which is give a global overview (but the information is heavily redundant). I could just pick one of the 3 redundant ones and skip the others, but then I am worried my presentation would be biased.

Any advice on how to properly do gene set enrichment analysis to find broader terms to describe and understand what's going on? Or alternative ways that lead to the goal of broader description and understanding? I would be very grateful.

Thank you!

# Answer
> 4 votes
Additionally there are plenty of tools to summarize long lists of GO terms by removing redundancy (which seems to be your concern), such as REVIGO (web) or rrvgo (R+Bioconductor). Disclaimer: I'm the author of the latter.

Such tools objectively pick the best representative term of each set of redundant terms, eg. by considering the GO term size (in genes), observed/expected ratio, p-value, etc.

# Answer
> 3 votes
It is totally acceptable in your main figure to display only a selected part of your GO term, as long as the reader is aware of this action (e.g., by including a legend specifying "selected pathways enriched"). Additionally, ensure that either readers have sufficient information to replicate your analyses or provide the full result table as supplementary information. By adhering to these points, you are transparent on what you did, allowing others to replicate the analyses independently. Purposely excluding top hit pathways that contradict your theory could be considered misleading.

If you are looking for broader terms, have you tried to a use slim GO? Those are GO with only a subset of the term. You could also look at the GO graph for your term and select the one higher in the graph. There are also plenty of tools out there to display GO results-graph. GOplot, GOnet are the first hit I got on google.

---
Tags: scrnaseq, single-cell, cancer, gsea
---


# Question
Title: Hybrid assembly versus polishing for hifi and illumina reads
I will have to carry out a project of assembly using hifi reads for which I have already illumina reads and I am wondering which of the hybrid assembly or polishing would be the best option for this kind of data? I wanted to use hifiasm which seems the gold standard and then nextpolish2 but I am wondering if using for example spades to do an hybrid assembly would lead to something better or worse? This question is only intended to see if someone did such benchmarking before or not.

Best regards,

# Answer
Usually, if you have HiFi reads, there is no way you can improve the assembly with short reads. Not just that they don't have higher quality, but they also map a lot worse. People sometimes use short reads for genome profiling and QCing the whole thing using Merqury.

Not sure if there is benchmarking published though.
> 1 votes

---
Tags: assembly, illumina, pacbio
---


# Question
Title: Why is my phylogenetic tree always splitting into one leaf and one node?
I used RAxML software on about 300 species mtDNA region (COI gene).
Unfortunately I used a highly variable region so my tree only has highly support values for branches close to the leaves.

But I also noticed that all the nodes split into one leaf and one other node... Why is this happening?
Isn't it common that a node will simply split into 2 other nodes?

I am using the `GTRGAMMA` function in RAxML (I don't know if that may be the cause).

Do you have any suggestions on this topic? Should I be suspicious of my tree if I don't see any node splitting into 2 other nodes?

**More info**

* Alignment algorithm: MUSCLE
* Graphical example:
My tree looks like the one on the left. Every bifurcation links to a tip/leaf (terminal node) of the tree and a non-terminal node of the tree. It sounds suspicious to me that I never observe a bifurcation that splits into 2 non-terminal nodes of the tree (figure on the right)

# Answer
What I suggest is the description sounds like a polytomy (star tree). In a phylogenetic tree, a polytomy is a node that has more than two immediate descending branches.

A polytomy can have many branches originate from one node. In contrast, a bifurcating tree only two branches stem from one node.

A polygomy simply means that the evolutionary relationships can not be fully resolved to dichotomies, i.e. bifurcations. There's nothing wrong with this result and it makes reading the tree much easier. The model you describe looks fine.

The only thing to watch if you're not looking at a dendrogram, but are looking at the phylogram. I doubt thats the explanation in any case.

---

From the update. The figure you presented is an asymmetry phylogeny and is perfectly normal, for example this is observed in molluscs, described in Palmer 1995. If a species is moving directionally e.g. east to west across a massive landmass and at uniform geographic intervals establishing an ecological niche and forming a new species this will result in an asymmetric tree and might form a "cline".
> 1 votes

---
Tags: phylogenetics, phylogeny
---


# Question
Title: Best reference book for basic concepts when working with bioinformatics
I'm starting to study Bioinformatics, I already have basic knowledge of Biochemistry and Genetics, but I want to buy a reference book, to use when I forget a concept (I'm a Computer Science student, not a Biotechnologist).

I have already had contact with the following two books (for basic concepts):

* Molecular Biology of the Cell, Alberts;
* Principles of Biochemistry, Lehninger.

Which is better for my use or do you have another recommendation?

I prefer a generalist book, which covers the basic concepts, and not a Bioinformatics book, which covers basic concepts in one chapter.

# Answer
For a reputable up to date text book on biochemistry I would only consider

> 'Biochemistry" Tenth Edition ©2023 by Jeremy Berg; Gregory Gatto Jr.; Justin Hines; John L. Tymoczko; Lubert Stryer published by MacMillan.

Stryer has been in the text book biochemistry business for a long time. I agree it's not a cheap option and I stopped reading biochemistry text books some time ago (see below), but this should be as good as it gets.

Don't use Lehninger although admittedly the 7th Edition is 2017, but biochemistry is a very fast moving field and you need to be as up to date as possible if you are relying on text books as a primary source of information. To be relying on text books they need to have a "Grade A" reputation and it to be extremely to date. Few text books stand the test of time.

I've taught text book stuff from courses that relied on "out of date editions", i.e. the authors have not updated the edition for a while. It can be a nightmare because essentially you need to navigate material that you know is incorrect and insubstantial, then accommodate and moderate these deficiencies in the end of course exam. Thats okay at a certain level, at research grade it gets unstuck.
> 0 votes

---
Tags: gene, computational-biochemistry
---


title: missing additional visits in GSE data

text:
# Question
Title: missing additional visits in GSE data
When downloading GSE data, often I encounter that the authors refer to additional visits in the data description, yet, they are not found in the data.

For example, the following study. They mention additional visits, however inspecting the phenotype data extracted using `getGEO` function does not show them (only baseline). The same happens with different studies as well, for example this study.

Example to how I extract the phenotype data:

```
gsm <- getGEO('GSE88885', GSEMatrix = TRUE)
gsm <- gsm[[1]]
pheno_dat <- gsm@phenoData@data %>% as.data.frame()

```

Am I missing some additional data?

# Answer
If it's not in there then you can either email the authors or just live with it. It is fully on the authors what you upload to NCBI. Nobody at NCBI checks the actual content of the uploads in terms of whether it indeed contains what methods texts or the paper pretend that it contains. Nor do (usually) the reviewers. Published data are great in terms of having a large repository of resources and knowledge available, downside is that metadata are often a mess or absent. There is no general and easy answer to this.
> 1 votes

---
Tags: rna-seq, differential-expression, public-databases, geoquery, gse
---


# Question
Title: How do I resolve disagreements in the determination of significant genes?
I want to find some good predictors (genes). I have checked the Spearman correlation of the expression of each of 23 genes with dependent variables (responders Vs non-responders and I saw only the 5 genes identified initially showed significant correlation with the dependent variable.

This is my data, log transformed RNA-seq:

```
              TRG    CDK6 EGFR  KIF2C CDC20
    Sample 1  TRG12  11.39 10.62  9.75 10.34
    Sample 2  TRG12  10.16  8.63  8.68  9.08
    Sample 3  TRG12   9.29 10.24  9.89 10.11
    Sample 4  TRG45  11.53  9.22  9.35  9.13
    Sample 5  TRG45   8.35 10.62 10.25 10.01
    Sample 6  TRG45  11.71 10.43  8.87  9.44
    ...

```

23 differentially expressed genes between responder patients (TRG12) and non-responder patients (TRG45).

1- I tested each of 23 genes individually in this code and each of them gives p-value \< 0.05 remained as a good predictor; For example for CDK6 I have done

```
> glm=glm(TRG ~ CDK6, data = df, family = binomial(link = 'logit'))
>
> summary(glm)

Call:
glm(formula = TRG ~ CDK6, family = binomial(link = "logit"),
    data = df)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.8286  -0.9317   0.6382   0.9500   1.5966

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  26.8950     9.6627   2.783  0.00538 **
CDK6         -2.6932     0.9765  -2.758  0.00582 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 77.347  on 55  degrees of freedom
Residual deviance: 64.062  on 54  degrees of freedom
AIC: 68.062

Number of Fisher Scoring iterations: 5

```

Finally I obtained five genes and I put them in this model:

```
final <- glm(TRG ~ CDK6 + CXCL8 + IL6 + ISG15 + PTGS2 , data = df, family = binomial(link = 'logit'))

```

but

```
> summary(final)

Call:
glm(formula = TRG ~ CDK6 + CXCL8 + IL6 + ISG15 + PTGS2, family = binomial(link = "logit"),
    data = df)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.5696  -0.9090   0.2233   0.7125   2.0131

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   3.2258    13.3877   0.241   0.8096
CDK6         -2.6101     1.0860  -2.403   0.0162 *
CXCL8         0.8702     0.9783   0.889   0.3737
IL6           0.9591     0.8883   1.080   0.2803
ISG15         0.9574     0.5101   1.877   0.0605 .
PTGS2        -0.4623     1.1005  -0.420   0.6744
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 77.347  on 55  degrees of freedom
Residual deviance: 53.917  on 50  degrees of freedom
AIC: 65.917

Number of Fisher Scoring iterations: 6

```

I am just seeing CDK6 as significant. I done cross validation

```
> cv.glm(df, final, K=nrow(df))$delta
[1] 0.2130050 0.2125748

```

Likely the delta values are not greatly differ (Model is right ???)

Then I tried Sensitivity and Specificity

```
> threshold=0.5
> predicted_values<-ifelse(predict(final,type="response")>threshold,1,0)
> actual_values<-final$y
> conf_matrix<-table(predicted_values,actual_values)
> conf_matrix
                actual_values
predicted_values  0  1
               0 16  7
               1 10 23
> sensitivity(conf_matrix)
[1] 0.6153846
> specificity(conf_matrix)
[1] 0.7666667

```

2- I went this way for step-wise method

```
full.model <- glm(TRG ~., data = df, family = binomial)

library(MASS)
step.model <- full.model %>% stepAIC(trace = FALSE)
coef(step.model)

> coef(step.model)
(Intercept)        CDK6        CHGA        CSF3       CXCL8       ERBB2         IL6       ISG15       KRT14       KRT17
-23770.6261   -469.3660    376.3528   -238.5697   1404.5905    590.2586    789.5295    842.2058    423.0402   -651.3403
     MAGEA1         OSM       PTGS2
   496.6933   -637.0583   -500.1040

```

In addition to 5 genes another genes remains as predictors but likely none of them are significant by Pr(\>|z|) \> 0.05

```
> summary(step.model)

Call:
glm(formula = TRG ~ CDK6 + CHGA + CSF3 + CXCL8 + ERBB2 + IL6 +
    ISG15 + KRT14 + KRT17 + MAGEA1 + OSM + PTGS2, family = binomial,
    data = df)

Deviance Residuals:
       Min          1Q      Median          3Q         Max
-1.317e-04  -2.100e-08   2.100e-08   2.100e-08   1.396e-04

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)  -23770.6  2922188.9  -0.008    0.994
CDK6           -469.4    64920.6  -0.007    0.994
CHGA            376.4    46191.2   0.008    0.993
CSF3           -238.6    30857.1  -0.008    0.994
CXCL8          1404.6   172558.7   0.008    0.994
ERBB2           590.3    71989.5   0.008    0.993
IL6             789.5   100130.7   0.008    0.994
ISG15           842.2   102934.0   0.008    0.993
KRT14           423.0    52565.5   0.008    0.994
KRT17          -651.3    79489.6  -0.008    0.993
MAGEA1          496.7    60319.6   0.008    0.993
OSM            -637.1    82121.8  -0.008    0.994
PTGS2          -500.1    73541.7  -0.007    0.995

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 7.7347e+01  on 55  degrees of freedom
Residual deviance: 1.4741e-07  on 43  degrees of freedom
AIC: 26

Number of Fisher Scoring iterations: 25

```

3- I went to lasso regression

```
> library(glmnet)
> lassoModel <- glmnet(
+   x=data.matrix(df[,-1]),
+   y=df$TRG,
+   standardize=TRUE,
+   alpha=1.0,
+   family="multinomial")
> #Derive coefficients for each gene to each subtype
> co <- coef(lassoModel, s=idealLambda, exact=TRUE)
> co
$`TRG12`
24 x 1 sparse Matrix of class "dgCMatrix"
                       1
(Intercept)  4.441322326
ALB         -0.010867331
AQP9        -0.051908908
CALML5       0.126359364
CCL4        -0.059869284
CDK6         0.283781437
CHGA        -0.088283426
CREB3L3      0.109988758
CSF3        -0.012311176
CXCL5       -0.030532485
CXCL6       -0.063610112
CXCL8       -0.084385732
ERBB2       -0.102303095
FGFRL1       0.148759702
IL1B        -0.072225606
IL6         -0.122671026
ISG15       -0.297386587
KLK5         0.039891341
KRT14       -0.004536929
KRT17        0.041693665
MAGEA1      -0.154075458
OSM          0.028987648
PTGS2       -0.034960405
S100A7A     -0.008190105

$TRG45
24 x 1 sparse Matrix of class "dgCMatrix"
                       1
(Intercept) -4.441322326
ALB          0.010867331
AQP9         0.051908908
CALML5      -0.126359364
CCL4         0.059869284
CDK6        -0.283781437
CHGA         0.088283426
CREB3L3     -0.109988758
CSF3         0.012311176
CXCL5        0.030532485
CXCL6        0.063610112
CXCL8        0.084385732
ERBB2        0.102303095
FGFRL1      -0.148759702
IL1B         0.072225606
IL6          0.122671026
ISG15        0.297386587
KLK5        -0.039891341
KRT14        0.004536929
KRT17       -0.041693665
MAGEA1       0.154075458
OSM         -0.028987648
PTGS2        0.034960405
S100A7A      0.008190105

>

```

Likely none of genes is significant so I just tried 5 genes from first glm

```
> finalLasso <- glm(df$TRG ~ CDK6 + IL6 + ISG15 + CXCL8 + CALML5 , data=df, family=binomial(link="logit"))
> summary(finalLasso)

Call:
glm(formula = df$TRG ~ CDK6 + IL6 + ISG15 + CXCL8 + CALML5, family = binomial(link = "logit"),
    data = df)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.6780  -0.6733   0.1643   0.6925   1.9182

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.4178    14.7456  -0.028   0.9774
CDK6         -2.2672     1.1537  -1.965   0.0494 *
IL6           1.2524     0.9197   1.362   0.1733
ISG15         1.1763     0.5479   2.147   0.0318 *
CXCL8         0.5413     0.6443   0.840   0.4008
CALML5       -0.7890     0.4219  -1.870   0.0615 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 77.347  on 55  degrees of freedom
Residual deviance: 49.839  on 50  degrees of freedom
AIC: 61.839

Number of Fisher Scoring iterations: 6

>

```

That 2 genes are significant; really confusing :(

Please somebody save me from this horrible confusion, finally which genes are significant remaining in model and which not

> First confusion: Why when I am putting each of 23 genes individually in glm model, 5 of them returns p-value \< 0.05 but when summary of final model containing these 5 genes says only p-value of one of them is \< 0.05
>
> Second confusion: When I am using stepwise regression, summary of model says none of my genes is significant while some of them had been remained as predictors by `coef(step.model)`
>
> Third confusion: CDK6 as my significant genes in two of model shows negative coefficient while I know that show be correlated with dependent variable positively

# Answer
confusion: Why when I am putting each of 23 genes individually in glm model, 5 of them returns p-value \< 0.05 but when summary of final model containing these 5 genes says only p-value of one of them is \< 0.05.

Ans. Pr(\>|z|) simply suggests that whether z-value deviates beyond average z-value and if so,what is the extent of its probability. less than .05 indicates a 95% confidence in computed z value i.e. an acceptable significance level. The summary of 5 genes simply tells significance of over-all model.
> 0 votes

---
Tags: r, rna-seq, linear-regression
---


# Question
Title: How to design a synthetic wastewater FASTQ file that could be flagged as "engineered"?
I am an apprentice-level python user who occasionally works in medical laboratories and mainly does business work. I am presented with the unenviable task of orchestrating a proficiency test of some sequencing teams who have the almost equally unenviable task of evaluating metagenomic FASTQ files of wastewater (sewage) data where the question is "could this have been engineered ?"

Trick: how do you make some FASTQ files that look like wastewater (plenty of those online), definitely has some pathogen (SARS-CoV-2, Ebola, whatever), oh, and the pathogen has to look like it could have been engineered? And make an equal number that do not have the pathogen but otherwise look the same?

Here is my initial thoughts:

1. World's simplest version: cat a virus FASTQ onto the end of a wastewater FASTQ.
2. change all the identifiers in the virus file so they look like a continuation of the wastewater, chop the virus file with visual mode in Vim so the reads are the same length as the wastewater.
3. instead of cat, do some python to build a new FASTQ by randomly selecting lines from the virus and wastewater "feedstock" FASTQs.
4. Very fancy: grab a assembled virus FASTA, splice in an interesting bit of pangolin, raccoon, pig, etc, then randomly select and subtly mutate substrings as reads of appropriate length, and assemble a de novo synthetic virus FASTQ. Mix that into the wastewater FASTQ as above.

Surely I am missing many, many steps and "very fancy" is probably not that fancy (or maybe the sun will grow cold before it succeeds, I do not know). Help and suggestions welcome, it could be fun.

# Answer
Use existing software to generate simulated fastq files. The input for these tools is generally a set of reference organisms that you'd expect to find in your sample. Simply add your 'engineered' viral genome to the list of references. These tools have error models to simulate all the errors that real sequencing instruments make.

How you'd go about making a viral sequence recognisable as genetically engineered, I have no clue. I'm sure the person who gave you this assignment has an idea on how to engineer viruses and what recognisable features that would result in, maybe consult with them.
> 1 votes

---
Tags: python, fastq, metagenome, quality-control
---


# Question
Title: How are BLAST record scores calculated in biopython
Can someone explain how the score is calculated for Bio.Blast.Record.Description? For example, if `Record.Description.num_alignments` is \>1, is `Record.Description.score` simply the sum of the alignment scores for the description's hit? Can alignments overlap? How does the score of a given alignment deal with the scores of overlapping HSPs in that alignment?

I'm getting the `Blast.Record` object via NCBIXML.parse(), and would like to rely on the description data as much as possible to remove the need for code that parses the Record.Alignments objects, but I'm a nervous including things when I don't know exactly how they were created.

Cheers! Tim

# Answer
> 1 votes
The results you'd get are the same as with any other NCBI BLAST+ tools, python does not compute it.

Biopython can run a blast query `qblast` on the internet or locally. To make it run locally you need to install blast tools on your machine and biopython will launch it as a subprocess. Here you can find the pertinent information: https://biopython.org/docs/dev/Tutorial/chapter\_blast.html#running-blast-locally

Filtering data is as much a burden as it is important. Depending on your question you want to answer. Some people prefer a low e-value. The score is however affected by the length of the sequence since the number of matches is bigger, more details here: https://www.nlm.nih.gov/ncbi/workshops/2023-08\_BLAST\_evol/blast\_score.html

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


# Question
Title: Abnormal representation of differentially expressed genes in volcano plots located at the extremities
*This question was also asked on reddit*

Does anyone know how to correct the abnormal genes that seemed to be located on the curves near the extremities? The analysis was done using R package edgeR on samples measured once. However, we couldn't figure out what was causing the anomaly or how to fix it. Does anyone have any advice on this matter? The code used for analysis is from the following website: https://blog.csdn.net/u012110870/article/details/102804557

# Answer
> 1 votes
In "normal" (that is replicated) experiments such structures usually are due to poor/insufficient/no prefiltering when there are unreliable count outliers (like two samples have zero counts for a gene and one has a lot).

As you cannot really meaningfully do any prefiltering other than requiring that either one or both samples have a certain read count you should probably just see these "DE results" (it's not differential really as it is unreplicated) as a ranking. Fold changes and p-values are very unreliable at 1 vs 1 so just rank by something like signed p-value and go along with that. Don't overinterpret these sorts of plots which are mostly meaningless in such a setup.

That having said, for prefilter I would probably require that at least one sample has a good number of reads, like 50 to at least dampen unreliably large fold changes a little bit in this setup.

# Answer
> 0 votes
Might as well transfer my reddit answer here...

If your p-value is smaller than the reciprocal of the number of atoms in the universe, the model used for calculating that differential expression is almost certainly wrong.

Do an MA plot; that will give you a better idea of what's going on.

* X axis - mean of (log expression for all samples); should be in the range of 0~15
* Y axis - log2 fold change
* colour - p-value (linear scale)

I suspect that the curves you're seeing are likely to be low count genes, in other words, genes on the left hand side of the graph in an MA plot.

I also suspect, given the super-astronomical p-values, that you've got an offset L2FC (i.e. the "spine" of the MA plot doesn't sit at L2FC = 0), suggesting that a covariate hasn't been appropriately accounted for. This is common when processing data from multiple samples from multiple different methods.

---
Tags: r, rna-seq
---


# Question
Title: How to download the same 18S rRNA gene region for multiple species?
I need 18S rRNA gene sequences for a wide variety of eukaryotic species.
I know that Ensembl has this gene sequence, but I don't understand what is the correct protocol to download the sequences. I already know the basics of Ensembl's Rest API, but there are still a few things that I don't understand.

I only know the gene name (18S rRNA) and the list of taxa I want to select it from. But I do not know the region of the gene.

For example: I tried looking for '18srrna' on the Ensembl browser and I noticed that for the Fruit fly there are 2 sequences corresponding to this gene (Ensembl IDs: FBgn0085813 and FBgn0267521).

1. How do I know which sequence I should use to make an alignment with other species?
2. How can I batch download the same region for all the species of interest?
3. Do you know any protocol/guide that goes into detail of how to download and select sequences to make phylogenies?

I don't have a solution for question (1) but for question (2) I am planning on requesting from the API the output of the endpoint `GET sequence/region/:species/:region`. But then again the problem is that I would need the position of the region for each species, and it sounds unfeasible to do for more than a dozen of species.

# Answer
> 0 votes
* How do I know which sequence I should use to make an alignment with other species?

Use both. For fruit flies there should thousands of 18S sequences.

* How can I batch download the same region for all the species of interest?

Using Blast (below)

* Do you know any protocol/guide that goes into detail of how to download and select sequences to make phylogenies?

Sure - Blast ... align via muscle -super5 and build via IQ-TREE

EMBL I don't know - using REST API via a CLI is possible. The problem is that you could potentially 1000s of 18S genes but a simple search term can ignore everything that isn't part of a genome explicitly labelled "18S". Thats where a Blast style approach is useful.

However, for NCBI Blast they intentionally block access to CLI with \>100 hits requested. This is a complete pain and really, really limits using Blast to extract out a gene remotely using a remote request script.

Basically the way to do it is to ditch the API and run this via the Blast NCBI web interface, or download the entire eukaryotic database and perform local Blast, MMseqs2 or Minimap2.

It's a bit bizarre I know, NCBI just doesn't like bots running its Blast service, but it gives more permission to use its web-service. The central advantage is that Blast will nicely remove the 18S gene from any given genome which can be downloaded as a single locus, without needing to download loads of genome stuff thats not needed. NCBI still does impose a limit on how much can be blasted via the web which sadly now is only 5000 (I'm sure it was much more like 50000) and its really not cool.

What you can do is use the Blast search "exclude/include" parameters to target the taxa of interest within the 5000 limit and using a tree of relevant eukaryotic evolution navigate the data sets. I assume you are looking for around 100000 sequences. Basically, once you've blasted you have the option of downloading just the region blasted and nothing else ... or the whole genome (avoid this option)... you will then have 5000 nicely 18S genes excised *in silico*

You can also set a lot other blast parameters, in particular only considering full length hits.

# Answer
> 0 votes
I can suggest two ways to batch-download aligned sequences for the species of interest. Both ways are integral to the Ensembl database and toolsets.

But before that, I suggest FBgn0085802 (which is rRNA) instead of FBgn0085813 and FBgn0267521 (which are pseudogenes). You may check FlyBase database.

1. Since you are interested in Fruit fly 18S rRNA homologs, Metazoa Ensembl is preferable. You query Ensembl ID in the Ensembl search box (at www.metazoa.ensebml.org) and click the popping link. 18SrRNA:CR41548 comes for FBgn0085802. On the left-hand side, Gene-based displays, select Sequence. under the Marked-up sequence, pick "BLAST this sequence".

Opening the BLAST/BLAT search automatically sets a DNA alignment query. You may search against the species of interest and add "18S" as Description. Later, run without changing any other parameter. It may take several minutes. Then, view results for hits found. You can investigate Sequence/s with the best alignment and closest genomic locations, then export data in FASTA or other formats.

While this assures you of the alignment and lets you investigate in detail the homologs, you may not batch-download particularly in this way.

2. As mentioned, BioMart Restful API is a way to batch-download sequences that you want. The website does only allow one search at a moment but editing the Perl script, which is generated according to your query, you may get sequences for multiple datasets. For setup, you should follow BioMart Perl API guideline cautiosly.

For your example query, set proper datasets, filters and attributes. Again for fruit fly, you may refer to:

https://www.metazoa.ensembl.org/info/data/biomart/biomart\_perl\_api.html

My example query seems like below in popping Perl script:

```
$query->setDataset("dyakuba_eg_gene");
$query->addFilter("with_dmelanogaster_eg_homolog", ["Only"]);
$query->addFilter("biotype", ["rRNA"]);
$query->addFilter("transcript_biotype", ["rRNA"]);
$query->addAttribute("ensembl_gene_id");
$query->addAttribute("ensembl_transcript_id");
$query->addAttribute("gene_biotype");
$query->addAttribute("description");
$query->addAttribute("transcript_biotype");

```

---
Tags: sequence-alignment, ensembl, multiple-sequence-alignment, dna
---


title: WGCNA soft threshold selection

text:
# Question
Title: WGCNA soft threshold selection
Here I am conducting WGCNA analysis. When calculating the soft threshold, it is basically below 0.8. I choose Power=6, which only has 0 and 7. Is it feasible?

# Answer
I don't know the context (how many samples, what tissue, experimental design etc) but the mean and median connectivities are pretty high at power 6. I would first check the data for strong drivers of correlation and possibly adjust for them if they do not seem biologically relevant to the question you study. Assuming you find no issues, I would go with a higher power here, maybe 12.
> 0 votes

---
Tags: wgcna
---


# Question
Title: How to add corresponding "chr" before location (start, end) of bed alike file?
I have some bed alike files but don't have `chr` term in the fisrt column, that I want to add `chr` and its corresponding chromosome number, I'm wondering how can I do that in bash!

```
head dat.tsv
9901    10000   0.186130156985236
10001   10100   244.668091357092
10101   10200   367.420929888855
10201   10300   179.801731647738
10301   10400   141.551984387272
10401   10500   92.5066880216621
10501   10600   92.3205578646769
10601   10700   127.96448292735
10701   10800   74.6381929510795
10801   10900   62.6327978255318

```

# Answer
> 1 votes
You could use the following perl one-liner, however if genome you are using contains more than one chromosome you need to figure out which lines get which chromosome chr1, chr2. If it's single chrosome like E.coli, usually name of chromosome is "chr".

```
perl -npe "s/(.*)/chr\t\1/" dat.tsv

```

# Answer
> -1 votes
I appreciate your responses.

Using the R `dplyr` package after using `bedtools`, I was able to deal with this problem. Using the `makewindows` function from `bedtool`, I first obtain a 100 bp window.

The `full_join` function from the `dplyr` R package then allows me to obtain a satisfactory result.

I was looking for a script that only uses bash.

---
Tags: bash, bed, tsv
---


# Question
Title: How can I check if a PDB file is incomplete?
One of the requirements of simulating the protein chain is to check for missing atoms or residues. I.e., using software tools to identify (and, if necessary, model) any missing regions.

How can I check if a PDB file is complete or incomplete using Python?

# Answer
> 3 votes
OK, lets say you need that information for `2PE4.pdb`, if you look at its webpage on RCSB PDB repository: 2PE4 and display its .pdb format file at https://files.rcsb.org/view/2PE4.pdb, you'll notice at `REMARK 465`:

```
REMARK 465
REMARK 465 MISSING RESIDUES
REMARK 465 THE FOLLOWING RESIDUES WERE NOT LOCATED IN THE
REMARK 465 EXPERIMENT. (M=MODEL NUMBER; RES=RESIDUE NAME; C=CHAIN
REMARK 465 IDENTIFIER; SSSEQ=SEQUENCE NUMBER; I=INSERTION CODE.)
REMARK 465
REMARK 465   M RES C SSSEQI
REMARK 465     GLY A   437
REMARK 465     HIS A   438
REMARK 465     HIS A   439
REMARK 465     HIS A   440
REMARK 465     HIS A   441
REMARK 465     HIS A   442
REMARK 465     HIS A   443

```

this should be the missing residues from the full peptide sequence of the protein that is listed as SEQRES:

```
SEQRES   1 A  424  ARG SER PHE ARG GLY PRO LEU LEU PRO ASN ARG PRO PHE
SEQRES   2 A  424  THR THR VAL TRP ASN ALA ASN THR GLN TRP CYS LEU GLU
SEQRES   3 A  424  ARG HIS GLY VAL ASP VAL ASP VAL SER VAL PHE ASP VAL
SEQRES   4 A  424  VAL ALA ASN PRO GLY GLN THR PHE ARG GLY PRO ASP MET
SEQRES   5 A  424  THR ILE PHE TYR SER SER GLN LEU GLY THR TYR PRO TYR
SEQRES   6 A  424  TYR THR PRO THR GLY GLU PRO VAL PHE GLY GLY LEU PRO
SEQRES   7 A  424  GLN ASN ALA SER LEU ILE ALA HIS LEU ALA ARG THR PHE
SEQRES   8 A  424  GLN ASP ILE LEU ALA ALA ILE PRO ALA PRO ASP PHE SER
SEQRES   9 A  424  GLY LEU ALA VAL ILE ASP TRP GLU ALA TRP ARG PRO ARG
SEQRES  10 A  424  TRP ALA PHE ASN TRP ASP THR LYS ASP ILE TYR ARG GLN
SEQRES  11 A  424  ARG SER ARG ALA LEU VAL GLN ALA GLN HIS PRO ASP TRP
SEQRES  12 A  424  PRO ALA PRO GLN VAL GLU ALA VAL ALA GLN ASP GLN PHE
SEQRES  13 A  424  GLN GLY ALA ALA ARG ALA TRP MET ALA GLY THR LEU GLN
SEQRES  14 A  424  LEU GLY ARG ALA LEU ARG PRO ARG GLY LEU TRP GLY PHE
SEQRES  15 A  424  TYR GLY PHE PRO ASP CYS TYR ASN TYR ASP PHE LEU SER
SEQRES  16 A  424  PRO ASN TYR THR GLY GLN CYS PRO SER GLY ILE ARG ALA
SEQRES  17 A  424  GLN ASN ASP GLN LEU GLY TRP LEU TRP GLY GLN SER ARG
SEQRES  18 A  424  ALA LEU TYR PRO SER ILE TYR MET PRO ALA VAL LEU GLU
SEQRES  19 A  424  GLY THR GLY LYS SER GLN MET TYR VAL GLN HIS ARG VAL
SEQRES  20 A  424  ALA GLU ALA PHE ARG VAL ALA VAL ALA ALA GLY ASP PRO
SEQRES  21 A  424  ASN LEU PRO VAL LEU PRO TYR VAL GLN ILE PHE TYR ASP
SEQRES  22 A  424  THR THR ASN HIS PHE LEU PRO LEU ASP GLU LEU GLU HIS
SEQRES  23 A  424  SER LEU GLY GLU SER ALA ALA GLN GLY ALA ALA GLY VAL
SEQRES  24 A  424  VAL LEU TRP VAL SER TRP GLU ASN THR ARG THR LYS GLU
SEQRES  25 A  424  SER CYS GLN ALA ILE LYS GLU TYR MET ASP THR THR LEU
SEQRES  26 A  424  GLY PRO PHE ILE LEU ASN VAL THR SER GLY ALA LEU LEU
SEQRES  27 A  424  CYS SER GLN ALA LEU CYS SER GLY HIS GLY ARG CYS VAL
SEQRES  28 A  424  ARG ARG THR SER HIS PRO LYS ALA LEU LEU LEU LEU ASN
SEQRES  29 A  424  PRO ALA SER PHE SER ILE GLN LEU THR PRO GLY GLY GLY
SEQRES  30 A  424  PRO LEU SER LEU ARG GLY ALA LEU SER LEU GLU ASP GLN
SEQRES  31 A  424  ALA GLN MET ALA VAL GLU PHE LYS CYS ARG CYS TYR PRO
SEQRES  32 A  424  GLY TRP GLN ALA PRO TRP CYS GLU ARG LYS SER MET TRP
SEQRES  33 A  424  THR GLY HIS HIS HIS HIS HIS HIS

```

Pay attention that the .pdb starts at residue position 20; that is a cloning artifact, as shown by:

```
SEQADV 2PE4 ARG A   20  UNP  Q12794              CLONING ARTIFACT
SEQADV 2PE4 SER A   21  UNP  Q12794              CLONING ARTIFACT
SEQADV 2PE4 THR A  436  UNP  Q12794              CLONING ARTIFACT
SEQADV 2PE4 GLY A  437  UNP  Q12794              CLONING ARTIFACT
SEQADV 2PE4 HIS A  438  UNP  Q12794              EXPRESSION TAG
SEQADV 2PE4 HIS A  439  UNP  Q12794              EXPRESSION TAG
SEQADV 2PE4 HIS A  440  UNP  Q12794              EXPRESSION TAG
SEQADV 2PE4 HIS A  441  UNP  Q12794              EXPRESSION TAG
SEQADV 2PE4 HIS A  442  UNP  Q12794              EXPRESSION TAG
SEQADV 2PE4 HIS A  443  UNP  Q12794              EXPRESSION TAG

```

See PDB residue numbering

Using **Biopython** lib you could try using `2pe4.cif` instead of .pdb file (https://files.rcsb.org/view/2PE4.cif) as input :

```
# First import MMCIF2Dict function
from Bio.PDB.MMCIF2Dict import MMCIF2Dict

# read in the CIF file.
pdb_info = MMCIF2Dict("2pe4.cif")

#
print(pdb_info["_pdbx_unobs_or_zero_occ_residues.label_seq_id"])
print(pdb_info["_pdbx_unobs_or_zero_occ_residues.label_comp_id"])
print(list(zip(pdb_info["_pdbx_unobs_or_zero_occ_residues.label_seq_id"],
pdb_info["_pdbx_unobs_or_zero_occ_residues.label_comp_id"])))

```

Result:

```
['418', '419', '420', '421', '422', '423', '424']
['GLY', 'HIS', 'HIS', 'HIS', 'HIS', 'HIS', 'HIS']
[('418', 'GLY'), ('419', 'HIS'), ('420', 'HIS'), ('421', 'HIS'), ('422', 'HIS'), ('423', 'HIS'), ('424', 'HIS')]

```

but for this kind of stuff it's probably better suited

PDBx Python Parser

---
Tags: python, pdb
---