cannin/biostars_qa huggingface dataset:
Question: **Question**:

What is the best and fastest way to query a VCF file for specific variants?

**Background**:

I have a tab-separated list of variants where the columns are CHROM, POS, REF, ALT. I would like to query a VCF file and get the records associated with only these specific variants. I know BCFtools, VCFtools, and tabix all allow you to supply a regions/positions file to search on CHROM and POS only, but I am interested in searching on CHROM, POS, REF, and ALT.

I also know this is very easy to do with grep, but grep doesn't take advantage of the VCF index file like the other tools do. As a result it is much slower, especially when searching very large VCF files.

Answer: Here's a solution that uses a combination of BCFtools, grep, and awk.

First, put your variants in tab-separated file (e.g. `variants.txt`) with columns `CHROM`, `POS`, `ALT`, and `REF`. For example:

1 215802298 A G
1 215844373 C T
1 215848808 A G
1 215901574 C T

Then query your VCF file (e.g. `my.vcf.gz`) like so:

variants='variants.txt'
vcf_in='my.vcf.gz'
vcf_out='variants_of_interest.vcf'

bcftools view -O v -R "$variants" "$vcf_in" \
| grep -Ef <(awk 'BEGIN{FS=OFS="\t";print "#"};{print "^"$1,$2,"[^\t]+",$3,$4"\t"}' "$variants") \
> "$vcf_out"

The output will be a VCF file with just the variants from your original list.

**Explanation:**

1. BCFtools quickly queries the VCF file based on `CHROM` and `POS` using the first two columns from `variants.txt`.
2. The `awk` command preprocesses `variants.txt` on-the-fly by adding necessary regular expressions:
- Add `#` to the variant list in order to preserve VCF header information in final results
- Prepend `^` to each variant in order to ensure that grep only searches from the start of each VCF line
- Add a placeholder regex for each variant's `ID` field
- Append a tab to each variant to ensure that grep does not include results for other indels in the VCF
3. Lastly grep filters the results from BCFtools based on `CHROM`, `POS`, `REF` and `ALT` to only include variants from the original variant list
Question: Hi,

I perform a single-cell analysis of mouse hematopoietic stem cells. I'd like to reveal the cell types within the dataset based on cluster-specific expression of marker genes, so I'm looking for a marker gene list of hematopoietic stem and progenitor cells.

More specifically, I'd like to identify the following cell types: LT-HSC, ST-HSC, MPP, LMPP, GMP, MEP and CLP.

Does anybody know where I can find such list for cluster annotation?

Thanks in advance.

Answer: [CellMarker][1] is a manually curated resource used for browsing, searching and downloading markers of diverse cell types of different tissues in human and mouse.

[1]: http://bio-bigdata.hrbmu.edu.cn/CellMarker
Question: Hi,

Anyone knows on which version of the human genome the data in the cbioportal is associated? I've downloaded mutations for several genes, but I've no idea whether they refer to hg18, hg19 or hg38, and didn't found any indication neither on the website nor on the forum.

thanks!

Answer: Sorry that this isn't explicit on the portal. All variants in all studies use `GRCh37.p13` loci and alleles. The autosomes and sex chromosome sequences are equivalent to `hg19`, but the mitochondrial DNA sequence is newer. [Click here](https://www.biostars.org/p/119295/#119308) for details.

The effect of each variant is mapped to canonical isoforms per gene ([as defined by Uniprot](http://www.uniprot.org/help/canonical_and_isoforms)) in [Gencode 19](https://www.gencodegenes.org/stats/archive.html#a19) using [Ensembl's VEP](http://useast.ensembl.org/info/docs/tools/vep/index.html) with the `vcf2maf` and `maf2maf` wrappers that [live here](https://github.com/mskcc/vcf2maf).
Question: When I run BaseRecalibrator, my input reference file is genome.fa, my mapped file is K1_bam.marked.realigned.fixed.bam and variants/dbSNP file is New.vcf. the commands are showed as follows:

java -jar GenomeAnalysisTK.jar -T BaseRecalibrator -I K1_bam.marked.realigned.fixed.bam -R genome.fa -knownSites New.vcf -o recal_dataK1.grp

Here is the result:

#####ERROR MESSAGE: Input files /raid1/rzeng/reference/New.vcf and reference have incompatible contigs: No overlapping contigs found.
#####ERROR /raid1/rzeng/reference/New.vcf contigs = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, X]
#####ERROR reference contigs = [chr10, chr11, chr12, chr13, chr14, chr15, chr16, chr17, chr18, chr19, chr1, chr2, chr3, chr4, chr5, chr6, chr7, chr8, chr9, chrM, chrX, chrY]

Here my reference contig order is obviously different from variants files. How can I adjust the order of reference file and make it the same as New.vcf file? Also, my .bam file is a mapped file by aligning sample sequencing reads with reference file, so .bam should has the same chromosome order as reference file, do I need to to adjust the chromosome order of .bam file too, So confused here!

Thank you

Answer: Easier is to sort your VCF instead of the BAM. As your BAM seems to be sorted on a simple lexical order you can just

sed 's/^\([^#]\)/chr\1/' < in.vcf | LC_ALL=C sort -t '\t' -k1,1 -k2,2n > out.vcf

beware if your VCF header contains some "##contig=" lines: you might have to change it too with sed .

update: added the sed

update **2021**: you can do this with `bcftools` only now.
Question: Hello,
I have 11 proteins in two datasets: with 5 in 1st and 6 in 2nd dataset. I have calculated the amino acid percentage for three amino acids (A, R, N) in them as shown in the attached file from excel ![enter image description here][1]

I want to normalize them by background amino acid probability of all proteins in two datasets and then conclude, for e.g., basic amino acids (R) are enriched or
hydrophobic amino acids (A) enriched. Kindly help.

Thank you!

[1]: http://i.imgur.com/X1pzo09.jpg

Answer: Percentages and frequencies are already a form of normalization, aren't they?

May be what you mean is how to test whether proportions of certain types of aa are significantly different in the two sets . A Fisher's exact test will do the job.
Question: I have a dds object in DEseq2 and I wanted to do row means filtering on my object. I have 30 total samples, therefore I have 30 total columns in the count table. I'd like to filter all read counts lower than 30 for each treatment (for rows 1-15 and rows 16-30). I've tried the following code:

keep = rowMeans(counts(dds[,c(1:15,16:30)])) >=30
dds_30=dds[keep,]

but I do not think this is working. Would anybody be able to help?

Answer: I'm not sure what your `dds` variable looks like, but this should work judging by your current code.

dds_30 <- dds[rowMeans(counts(dds)[, 1:15]) >= 30 & rowMeans(counts(dds)[, 16:30]) >= 30, ]
Question: Hi All,

Just a quick question re tophat- How does tophat handles the paired-end data when both the reads (left and right) are in one file. I mapped this data (still running) using tophat without separating the paired-end reads in two different files. Will top handle this or I should separate the files?

Thanks,RT

Answer: You should separate the FASTQ reads into two files for Tophat. It will not recognize them as paired-end reads if they are together in the same file.

From the manual:

```
##Using TopHat

Usage: `tophat [options]* [reads1_2,...readsN_2]`

When running TopHat with paired reads it is **critical** that the `*_1` files an the `*_2` files appear in separate comma-delimited lists, and that the order of the files in the two lists is the same.
```
Question: I have a slightly unusual sequencing set-up (from a Novaseq paired-end sequencing experiment). I have i7 barcodes which are the sample indices. My i5 index is an 9bp UMI of random N's. My goal is to demultiplex based on the i7 sample indices, and then extract a separate file that is i5 UMI fastq. I want to use umi-dedup or something similar to then annotate R1/R2 headers with the UMI.

I can extract the R1, R2, and I5 index files to a fastq using the following command and RunInfo.xml file

bcl2fastq --output-dir BCL --sample-sheet SampleSheet.csv --create-fastq-for-index-reads

but this gives me an index file with only N's

@A00454:609:H7VCGDRXY:1:2101:3821:1000 2:N:0:NTTACTCG
NNNNNNNN
+
########

Is there a way to somehow convert the N's to the actual UMI so that I can use it in downstream analysis? Any help would be appreciated since I'm very new to using UMIs in this way.

Answer: You need to add the option `--mask-short-adapter-reads 0` when you demultiplex the data with `bcl2fastq`. This will restore the sequence in your UMI file. Normally if a read falls below 22 bp it is masked by default, which is what you are observing.
Question: Hi,

I am using **samtools markdup -r** parameter to remove PCR duplicates from my mapped ChIP reads. How does it detects the PCR duplicated reads?

Please let me know if I am understanding it correctly. This is what I understood:

I need to name sort the mapped reads and put MS-MC tags (using fixmate) and then I need to coordinate sort the reads before using the "**markdup -r**" command.
If the chromosome number and 5' end of the read are the same, it will treat them as duplicate and by looking at the MC and MS scores it will choose the best-read among the duplicates and delete the rest of the reads (I am using **"-r"** parameter).. Is this what is happening here?

Any kind of help will be highly appreciated.

Thanks

-S

Answer: Samtools markdup is written to match Picard 2.10.3 (also Biobambam's bamstreamingmarkduplicates) so if you can find documentation on those then it should also apply to Samtools.

It may seem like a complex dance to have both name and position sorted requirements, but this is perhaps due to a traditional view and attempting to do markdup as an after thought.

The design leverages the fact that the data produced by the sequencing instrument and subsequent aligner will be name collated (not sorted, but collated is all it actually needs - we should probably fix the man page). So while it is in this order we can add additional tags regarding the other ends of reads in an efficient manner. Then after positional sorting markdup has sufficient information to make the assessment which it can now also do in an efficient manner.

When viewed in that perspective, marking duplicates actually has zero additional sort steps than we're already using in our standard pipelines. Hopefully that helps to explain the somewhat unusual approach.

As for the exact maths, I don't know all the intricacies, but I think it's pretty much as you describe. There is more information here:

https://www.htslib.org/doc/samtools-markdup.1.html

https://www.htslib.org/algorithms/duplicate.html
Question: I have the following FASTA file, `original.fasta`:

```
>foo
GCTCACACATAGTTGATGCAGATGTTGAATTCACTATGAGGTGGGAGGATGTAGGGCCA
```

I need to change the record id from `foo` to `bar`, so I wrote the following code:

```py
from Bio import SeqIO

original_file = r"path\to\original.fasta"
corrected_file = r"path\to\corrected.fasta"

with open(original_file) as original, open(corrected_file, 'w') as corrected:
records = SeqIO.parse(original_file, 'fasta')
for record in records:
print record.id # prints 'foo'
if record.id == 'foo':
record.id = 'bar'
print record.id # prints 'bar' as expected
SeqIO.write(record, corrected, 'fasta')
```

We printed the record id before and after the change, and get the expected result. We can even doublecheck by reading in the corrected file again with BioPython and printing out the record id:

```py
with open(corrected_file) as corrected:
for record in SeqIO.parse(corrected, 'fasta'):
print record.id # prints 'bar', as expected
```

However, if we open the corrected file in a text editor, we see that the record id is not `bar` but `bar foo`:

```
>bar foo
GCTCACACATAGTTGATGCAGATGTTGAATTCACTATGAGGTGGGAGGATGTAGGGCCA

```

We can confirm that this is what is written to the file if we read the file using plain Python:

```py
with open(corrected_file) as corrected:
print corrected.readlines()[0][1:] # prints 'bar foo'
```

Is this a bug in BioPython? And if not, what did I do wrong and how do I change the record id in a FASTA file using BioPython?

Answer: You have to change the record.description to achieve your goal:

```py
from Bio import SeqIO

original_file = "./original.fasta"
corrected_file = "./corrected.fasta"

with open(original_file) as original, open(corrected_file, 'w') as corrected:
records = SeqIO.parse(original_file, 'fasta')
for record in records:
print record.id
if record.id == 'foo':
record.id = 'bar'
record.description = 'bar' # <- Add this line
print record.id
SeqIO.write(record, corrected, 'fasta')
```
Question: Hi,

I have done survival analysis on a gene, using the coxph call and saved it in a variable a. I know i can do summary(a) and get a return with the coef, exp(coef), se(coef) and power- p-value. I want to take JUST the p-value out of a to put into a a matrix.

So essentially is there a way of accessing just the p-value form a coxph object , e.g. coefVal <- a$coef gives me the coef value.. ?

Thanks,
R

Answer: Assign coef summary to a variable

coeffs <- coef(summary(a))

coeffs should be something like this:

Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.86370076 0.71438758 -5.40840977 6.358681e-08
A 20.62019096 264.13705493 0.07806626 9.377753e-01
B 4.00226650 0.71574200 5.59177257 2.247631e-08
C 4.50171396 0.71849626 6.26546610 3.717119e-10
...

Subset 4th column as P is 4th column in this example:

as.matrix(coeffs[,4])
[,1]
(Intercept) 6.358681e-08
A 9.377753e-01
B 2.247631e-08
C 3.717119e-10
...
Question: Dear all,

I would like to write shell script to automatic alignment my pair end data with bowtie (or bwa) aligners.

Problem is, when I have lot of pair-end fastq files in one directory I have problem to loading data in for cycle.

Here is my code:

for i in $(ls *.fastq | rev | cut -c 22- | rev | uniq)
do ./bowtie2 --very-sensitive -p32 --rg-id ${i} --rg-sample ${i} --rg-library Illumina --rg-platform Illuimna -x $reference -1 ${i}_L001_R1_001.fatsq.gz -2 ${i}_L001_R2_001.fatsq.gz -S ${i%.fastq.gz}${i}.sam
done;

This is very nice solution help solved by RamRS. Is there any other way how to treat Bowtie or BWA with many samples (corresponding one sample to each read) in one for cycle?

my directory contains lot of fastq.gz:

name1_R1_001.fatsq.gz
name2_R2_001.fatsq.gz
...
...
...

Each sample name (prefix before `R._001.fastq.gz`) has different name length so algorithm above doesn't work properly :-(

Thank you for any idea and solutions..

Answer: ```
for I in `ls *.fastq | cut -d "_" -f 1` ;
do ./bowtie2 --very-sensitive -p 32 --rg-id $i --rg-sample $i --rg-library Illumina --rg-platform Illuimna -x $reference -1 $i\_L001_R1_001.fastq.gz -2 $i_L001_R2_001.fastq.gz -S $i.sam ;
done ;
```

edit: removed useless grep
Question: I have defined 44 linkage blocks as "interesting" through a GWAS - NGS cohort.

I have run certain tests on those, but I wish to understand how typical or atypical those results are by comparing the 44 loci to 44 similarly size loci that are "matched" according to various things.

1. Matched based on the number of genes and total gene content

2. Matched based on the number of regulatory elements (e.g. Conserved TFBS)

Is there a software that does this? If not, how to people do it?

Thank you.

Answer: Assuming your datasets are in sorted BED format, you could use `bedmap` with the `--count` operand to calculate counts of mapped genes and regulatory elements over a reference set (docs). Other statistical operations are available, score and non-score-based. You could use this to calculate statistics for treatment and control sets, observed and background sets, etc. Alternatively, if you don't already have a background set, you could, say, map over windows to find regions that have characteristics you want.
Question: Hi there I'm trying to follow the instructions here

[https://www.htslib.org/workflow/fastq.html][1]

However when I get to the second step and input:

samtools fixmate -O bam,level=1 -m my.sam fixmate.bam

It tells me that level is an unknown option.

I tried stripping that part of the command so it ran

samtools fixmate -O -m my.sam fixmate.bam

However it just tells me:

Usage: samtools fixmate
Options:
-r Remove unmapped reads and secondary alignments
-p Disable FR proper pair check
-c Add template cigar ct tag
--input-fmt-option OPT[=VAL]
Specify a single input file format option in the form
of OPTION or OPTION=VALUE
-O, --output-fmt FORMAT[,OPT[=VAL]]...
Specify output format (SAM, BAM, CRAM)
--output-fmt-option OPT[=VAL]
Specify a single output file format option in the form
of OPTION or OPTION=VALUE
--reference FILE
Reference sequence FASTA FILE [null]

As elsewhere in samtools, use '-' as the filename for stdin/stdout. The input
file must be grouped by read name (e.g. sorted by name). Coordinated sorted
input is not accepted.

I tried transposing the -O and -m switches out of desperation but that doesn't run either.

Does anybody perhaps know what I'm doing wrong? Thanks

[1]: https://www.htslib.org/workflow/fastq.html

Answer: So my problem was that I'd installed samtools from the universe repository and it was a really old version.

Now I installed the latest version manually the instructions I was following actually work fine.

Damn you universe repository guy!
Question: Hello everyone,
I am new to perl API for EnsEMBL. I am even new to perl but doing some readings on both of these, I could write a small code to extract the protein and cDNA sequences of canonical transcript for given Ensembl gene ids.
I have pasted my code below:

#!/usr/bin/env perl

use strict;
use warnings;

use Bio::EnsEMBL::Registry;
use Data::Dumper;

my $registry = "Bio::EnsEMBL::Registry";
$registry->load_registry_from_db( -host =>'ensembldb.ensembl.org', -user => 'anonymous' );
my $gene_adaptor = $registry->get_adaptor( 'Human', 'Core', 'Gene' );
my @gene_ids=("ENSG00000001629","ENSG00000001630","ENSG00000001631","ENSG00000002016");

foreach my $gene_id (@gene_ids) {

my $gene = $gene_adaptor->fetch_by_stable_id($gene_id);
my $cdsseq = $gene->canonical_transcript()->translateable_seq();
my $protseq = $gene->canonical_transcript()->translate()->seq();
print $gene_id,"_cDNA\n",$cdsseq,"\n",$gene_id,"_protein\n",$protseq,"\n";
}

But, what I really want to do is to get all the translatable transcripts (having protein sequence) for a gene including the ensembl protein id and transcript id. (ENSG..............|ENSP...................|ENST...............\n peptide sequence) in a fasta file. And similary, cDNA sequences in another fasta file. Also, my script breaks when there is no translatable sequence for a gene, how do i set an if clause for that?
Thank you for your help!!

Abhishek

Answer: Hello, you may probably check if the translatable sequence is define before your print. If the translatetable sequence is not define, you will not print values and go to next turn. Or if the translatable sequence is define, you print your values but do nothing if not.

For exemple:

#!/usr/bin/env perl

use strict;
use warnings;

use Bio::EnsEMBL::Registry;
use Data::Dumper;

my $registry = "Bio::EnsEMBL::Registry";
$registry->load_registry_from_db( -host =>'ensembldb.ensembl.org', -user => 'anonymous' );
my $gene_adaptor = $registry->get_adaptor( 'Human', 'Core', 'Gene' );
my @gene_ids=("ENSG00000001629","ENSG00000001630","ENSG00000001631","ENSG00000002016");

foreach my $gene_id (@gene_ids) {
if ( defined ( $gene->canonical_transcript()->translateable_seq() ) ) {
my $gene = $gene_adaptor->fetch_by_stable_id($gene_id);
my $cdsseq = $gene->canonical_transcript()->translateable_seq();
my $protseq = $gene->canonical_transcript()->translate()->seq();
print $gene_id,"_cDNA\n",$cdsseq,"\n",$gene_id,"_protein\n",$protseq,"\n";
}
}

I didn't test this short script but it may works.

Nolwenn
Question: Hello,

I am doing mouse RNA-Seq analysis. I want to visualize mapped bam file on IGV. I am having hard time downloading my bam files from linux remote machine to local computer. Can I visualize this directly on remote computer ?

I tried to read user guide from [IGV][1], but I still could not figure out how to run my commands ?

If you can recommend me any tutorials, I greatly appreciate it.

Thank you very much.

[1]: https://igv.org/

Answer: There are three ways to do this:

1. X forwarding
2. Remote Desktop
3. Mount the remote machine as a network share on the local machine.

In both the first two cases you need to install IGV of the remote linux machine. Older version are available from most of the main linux repos, so you should be able to install with something like `apt-get install igv` or `yum install igv`. Alternatively if you use the (ana)conda package management system on the remote machine, you can install a recent version with `conda install igv -c bioconda`.

# X forwarding

If you are using Windows or OSX on the local machine you will need to install an X server on your local machine. There are several windows X11 servers, that are either intergrated with terminal emulators or stand alone. The [VcXsrv][1] server is popular as a standalone solution and [mobaXterm][2] integrates a terminal and ssh solution.

In OSX you'll need to install [XQuartz][3].

Now connect to your remote server with your favourite SSH solution making sure you configure it to forward X11 connections. On the commnad line thats `ssh -X`. Start IGV on the remote machine and if your X Server is correctly configured, IGV will appear on your local machine.

# Remote Desktop
If the remote machine is a desktop linux machine, you should be able to run a remote desktop sharing daemon on it and then connect to that remote desktop from the local machine. I've successfully used [Chrome Remote Desktop][4] for this, but many other solutions exist.

# Mount the remote machine on the local machine

Basically in this solution you make the disk of the remote machine appear like a network drive on the local machine. The most fail safe way to do this is using `sshfs`, which use can use to connect to any machine that will accept incoming ssh connections.

[1]: https://sourceforge.net/projects/vcxsrv/
[2]: https://mobaxterm.mobatek.net/download.html
[3]: https://www.xquartz.org/
[4]: https://remotedesktop.google.com/headless/