cannin/biostars_qa huggingface dataset:
Question: Hello,

I have 15 vcf files for one type of population and 15 vcf files for another type. I want to check the differences between the two, and also the similarities. What changes from one group to another and what remains the same, and a signifcance score if possible.

I have read about PLINK but I am not sure how the pipeline should be. Which steps should I folllow? I read the documentation and it is not clear to me.

I also read about bcftools isec: which is useful to intersect multiple vcf files. So I could merge the 15 vcf files between them and the other 15 vcf files between them and end up with two files: population1_variants.vcf and population2_variants.vcf, and then compare those two against eachother and check for the differences and similarities?

Which approach is better? Is this the way people usually analyze variants among populations? How can I asess significance of the results? Are there any other approaches?

Thank you

Answer: It really depends on what you want to achieve with this comparison. You could merge all VCFs and do an association analysis between the two populations using plink to find differences between the two groups or you could do a PCA using all samples to see if the two populations have a clear separation between them.

Try doing an association analysis:

plink --file mydata --assoc

Look for SNPs with statistical significance between the two groups.

http://zzz.bwh.harvard.edu/plink/anal.shtml
Question: Through vcftools I got a file(`my.var-final.vcf`, 27 MB) which contain in formation of SNPs and indels:

```
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT my-sorted.bam
comp904_c0_seq1 30 . G T 73.5 . DP=4;VDB=0.0014;AF1=1;AC1=2;DP4=0,0,4,0;MQ=60;FQ=-39 GT:PL:GQ 1/1:106,
12,0:21
comp904_c0_seq1 37 . C T 52 . DP=4;VDB=0.0014;AF1=1;AC1=2;DP4=0,0,3,0;MQ=60;FQ=-36 GT:PL:GQ 1/1:84,9
,0:16
comp904_c0_seq1 41 . A T 64.3 . DP=6;VDB=0.0020;AF1=1;AC1=2;DP4=0,0,5,0;MQ=60;FQ=-42 GT:PL:GQ 1/1:97,1
5,0:27
comp904_c0_seq1 74 . A G 4.77 . DP=21;VDB=0.0147;AF1=0.4999;AC1=1;DP4=10,5,3,1;MQ=60;FQ=6.99;PV4=1,1.2e-06,1,1
GT:PL:GQ 0/1:33,0,255:33
comp904_c0_seq1 133 . G T 137 . DP=36;VDB=0.0404;AF1=0.5;AC1=1;DP4=2,3,19,10;MQ=60;FQ=33;PV4=0.35,1.6e-09,1,1
GT:PL:GQ 0/1:167,0,60:63
```

Is there any way to summarize this variation information, like some tools, scripts or something?

Answer: There are tons of tools that will give you what you want. Here is one: http://vcftools.sourceforge.net/documentation.html#file
Question: Hi,

I wish to perform variant analysis for exon capture data. The samples are human tissue samples. I already have the sequencing data and the list of exons to analyze. The exon co-ordinates are taken from the GCF_000001405.38 assembly.

I looked up on the UCSC genome browser and there the hg38 assembly version is GCA_000001405.15. So I looked up on NCBI and found this link to the latest version of the assembly.

[GRCh38.p12][1]

In this, I have downloaded the GCF_000001405.38_GRCh38.p12_genomic.fna.gz file. Is this the correct reference or should I use the one on UCSC genome browser?

[UCSC Genome Browser][2]

[1]: ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.38_GRCh38.p12/
[2]: http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz

I will appreciate the guidance.

Thanks!!!

Answer: GCA_000001405.15 and GCF_000001405.26 are identical. They are GRCh38 without any patches. Here's an explanation about patches from GRC: https://www.ncbi.nlm.nih.gov/grc/help/patches/
GCF_000001405.38 is the latest version of GRCh38 with patches. Depending on your use case, you may not need the assembly with the latest patches. In such a situation, you can use GRCh38 without any patches (GCF_000001405.26) from here: ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.26_GRCh38/GCF_000001405.26_GRCh38_genomic.fna.gz
Question: Hi,

I would like to download the FASTA files of all genes without an intron.

Answer: In R you could do something like this to get all transcript with only one exon (thus intronless)

library(biomaRt)
library(dplyr)

ensembl = useEnsembl(biomart="ensembl", dataset="hsapiens_gene_ensembl")
# will give you one line per exon
gene <- getBM(attributes=c('ensembl_gene_id','ensembl_transcript_id','hgnc_symbol','ensembl_exon_id'), mart = ensembl)

# extract all transcript with only one exon
gene.oneExon <-
gene %>%
group_by(ensembl_transcript_id) %>%
summarise(n=n()) %>%
filter(n==1)

# get sequence of associated transcript
gene.oneExon.seq <- getSequence(id=gene.oneExon$ensembl_transcript_id,type="ensembl_transcript_id", seqType="gene_exon",mart=ensembl)
Question: I'm just a beginner in this area so I apologize first for wasting your time with this stupid question:

I recently downloaded some data from Uniprot and want to do some domain coverage analysis, but just realized that there are some overlapping ranges that needs to be merged.

For example, my data looks like this:

no.	pfamseq_acc	seq_version	crc64	md5	pfamA_acc	seq_start	seq_end
1857719	I1HYZ2	1	E9A7BABADB78EF7A	217559545415239f73152e983a1e4226	PF00083	11	244
1857720	I1HYZ2	1	E9A7BABADB78EF7A	217559545415239f73152e983a1e4226	PF00083	482	737
1857722	A8I2W9	1	125E4A38BE3390D8	3be8943f5916ee0fd2c80adba0af3377	PF00083	83	299
1857723	A8I2W9	1	125E4A38BE3390D8	3be8943f5916ee0fd2c80adba0af3377	PF00083	286	500
1857724	A0A0L0MB57	1	277707E4DDF29006	cf70dbe1ad2c63df7e4b6d419c5dff05	PF00083	10	222
1857725	A0A0L0MB57	1	277707E4DDF29006	cf70dbe1ad2c63df7e4b6d419c5dff05	PF00083	223	423

And what I want is:

no.	pfamseq_acc	seq_version	crc64	md5	pfamA_acc	seq_start	seq_end
1857719	I1HYZ2	1	E9A7BABADB78EF7A	217559545415239f73152e983a1e4226	PF00083	11	244
1857720	I1HYZ2	1	E9A7BABADB78EF7A	217559545415239f73152e983a1e4226	PF00083	482	737
1857722	A8I2W9	1	125E4A38BE3390D8	3be8943f5916ee0fd2c80adba0af3377	PF00083	83	500
1857724	A0A0L0MB57	1	277707E4DDF29006	cf70dbe1ad2c63df7e4b6d419c5dff05	PF00083	10	423

I have tried to write it in R which looks like:

library(dplyr)
setwd("C:/Users/hy/Documents/Desktop/test2")
d<- read.csv("test2.tsv", header = T, sep = "\t", dec = ".")

d %>%
group_by(md5) %>%
arrange(md5,seq_start) %>%
mutate(indx = c(0, cumsum(as.numeric(lead(seq_start)) >
cummax(as.numeric(seq_end)))[-n()])) %>%
group_by(md5,pfamseq_acc,pfamA_acc)
summairse(start = min(seq_start), end = max(seq_end))

The "summarise" could give me the correct result, but I failed to write it in a new file.

So I would appreciate it if anyone could give some ideas about how to achieve that (preferably by R or Python), any other methods is also welcomed.

Thanks.

Answer: You could use `bedops --merge` to merge intervals in a sorted BED file:

$ bedops --merge in.bed > out.bed

A BED file is just a tab-delimited text file, minimally three columns, which specify chromosome, start, and stop attributes of each interval, respectively. You can export such a file from R via `write.table` or the like.

Using a toolkit like this will likely be faster than using a library in R or Python.
Question: Dear Friends,

I have a list of 8000 samples in a file "samples.txt":

samples.txt:
TCGA..barcode..
TCGA..barcode..
.
.

I am using bcftools to only keep these samples in the vcf.gz file. The vcf.gz file has 10000 samples. Hence, I am trying to use bcftools to keep only the 8000 samples in the "samples.txt" file in the vcf.gz file and remove the remaining 2000 samples. I did:

bcftools -S samples.txt vcf.gz -o filtered-vcf.vcf

it gives me error:

[E::main] unrecognized command -S

Could you please suggest me what could be the issue here, and how you think I can do the above? Thanks much.

Answer: subcommand 'view' is missing:

bcftools view -S samples.txt -o filtered-vcf.vcf vcf.gz
Question: I was looking into **ClinVar** data for getting mutation lists. There were mutations which were in the form `GENE:c.***` representing they are CDS mutations and `GENE:p.***` representing the amino acid changes.

**What are those in the following forms represent**?

1. `m.***`
2. `GENE:n.***`
3. `GENE:g.***`
4. `nsv***`

***Example***:

- `TBC1D24:c.1143-6C>T` - **CDS mutation**
- `NP_002760.1:p.Cys139Ser` - **Protein mutation**
- `m.1606G>A` - **??**
- `U43746.1:n.2241A>G` - **??**
- `NC_000023.11:g.53254331_53296102dup41772` - **??**
- `nsv513787` - **??**

Answer: They are the HGVS nomenclature.

The website is: http://www.hgvs.org/mutnomen/

But curiously, it is down right now.

[Here][1] is a cached copy.

[1]: http://webcache.googleusercontent.com/search?q=cache:4qlzlxEoo2kJ:www.hgvs.org/mutnomen/recs.html+&cd=1&hl=en&ct=clnk&gl=us
Question: I am working with [*rehh* package][1].

I have to create an object of class haplohh from the function data2haplohh:

test<-data2haplohh(hap_file="test.hap", map_file="test.map.inp", haplotype.in.columns=TRUE)
Map file seems OK: 1054416 SNPs declared for chromosome 1
Haplotype are in columns with no header

> Error in scan(hap_file, what = "character", quiet = TRUE) : too many
> items

So the *.hap file has 1054416 lines and 5008 columns(approx 10 GB size), *.map.inp file has 1054416 lines and 5 columns(15 MB size). My systems working RAM memory is 64 GB and 8TB is the storage space

I have an option of chopping up the files but, in this case It would create errors in my calculation. So I need to load the input files as a whole.

Please tell me how to solve this!

[1]: https://cran.r-project.org/package=rehh

Answer: Indices for R vectors are 32-bit signed integers, even in 64-bit versions of R, which can put an upper limit on how many elements you can work with.

Looking at the source of `getAnywhere(data2haplohh)` for your error message, it looks like it tries to make a vector or matrix with 5.3B elements (1054416 rows x 5008 columns). The limit you are perhaps running into is 2^31-1 or 2.1B elements.

Maybe contact the developers for help with adjustments to code to support larger inputs, or look at subsampling or validating your input.
Question: Are there any common approaches to make use of unmapped reads to reconstruct strain specific genome segments?

Are these just your normal denovo aproaches do you somehow also take the mapped reads and the reference in to account? Maybe for example read pairs in which one reads maps and not the other? Or reference assisted assembly?

And is this at all possible with Solid (fragment and PE ) sequencing data? A denovo assembly or reference assisted assembly?

Answer: There was a recent paper which did this with a ton of rice accessions:

[Exploring the rice dispensable genome using a metagenome-like assembly strategy][1]

>We proposed an assembly strategy incorporating SOAPdenovo [21], Meta-IDBA [22], and GICL [23] to assemble the dispensable sequences (“Materials and methods”; Fig. 1a). Both reads of a read pair of which at least one read could not be mapped to the reference genome were extracted to do de novo assembly (“Materials and methods”). In total, 395,891,829 read pairs (7.0 % of all the read pairs) were collected, of which 628,163,124 reads (79.3 %) could not be mapped to the Nipponbare genome

[1]: http://Exploring%20the%20rice%20dispensable%20genome%20using%20a%20metagenome-like%20assembly%20strategy
Question: Hi,

I'm trying to detect primer presence in my data via FastQC.

What's the direction of the graphs given by FastQC, 5'->3' or 3'->5'?

I expect my primers to be in the 3' end.

Many thanks

Answer: The plots are all 5' to 3'.
Question: I'm trying to automatically compare BAM files being output by bowtie2 (for a continuous integration system). Currently using `bamUtil diff` which works well, except I'm noticing sometimes 2 lines in the BAM file will have the same name/start position, but represent different reads. For whatever reason, bowtie2 sometimes reverses the order of the lines, so bamUtil's algorithm calls a mismatch. I see several potential solutions to this problem:

1. Get samtools to sort by another dimension in addition to position/name. Pretty much anything should do. Can this be done? I don't see any options in the man page.

2. Use another comparision tool than bamUtil which doesn't have this problem. Honestly a stand-alone program for diffing would be best for me anyway and I'm considering writing one.

3. Find some way to get bowtie2 to be more deterministic in its ordering. Any ideas here?

Answer: Exactly which parts of the alignments do you want to be diffed? If you just care about name/sequence then you can use bamHash. If you want more than that then you'll need a custom solution. There are no options to have samtools sort in a 3rd way (no one has ever needed that).
Question: Some of the SNPs in my dataset are without rs ids, so in the bim file, they have dots instead of rs ids for them. So for errorless analysis how should replace the rs id column with that of the chromosome: location column instead?

Answer: Using **awk**, (not tested):

awk '{if($2 == ".") {print $1,$1"_"$3,$4,$5,$6} else {print $0}}' myfile.bim > myfileClean.bim
Question: I went through the vignette about interaction terms and would like to understand if I am applying interaction terms correctly and biologically, if I am extracting the correct the terms. Here is the example code from DESeq2 R documentation for **two conditions (A, B) and three genotypes (I,II,III)**. Let A and B be 2 conditions as: **placebo and treatment** respectively and in case of genotype: **I is KO (where, gene X is knock out), II is WT and III is Tg (where, extra copy of gene X is present).**

> dds <- makeExampleDESeqDataSet(n=100,m=18)
> dds$genotype <- factor(rep(rep(c("I","II","III"),each=3),2))
> dds
class: DESeqDataSet
dim: 100 18
metadata(1): version
assays(1): counts
rownames(100): gene1 gene2 ... gene99 gene100
rowData names(3): trueIntercept trueBeta trueDisp
colnames(18): sample1 sample2 ... sample17 sample18
colData names(2): condition genotype
> dds$genotype <- factor(rep(rep(c("I","II","III"),each=3),2))
> design(dds) <- ~ genotype + condition + genotype:condition
> dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
> resultsNames(dds)
[1] "Intercept" "genotype_II_vs_I" "genotype_III_vs_I"
[4] "condition_B_vs_A" "genotypeII.conditionB" "genotypeIII.conditionB"

If we go through each and every terms under `resultsName(dds)` then:
> `genotype_II_vs_I` will be the comparison between WT and KO across treatment and placebo,

> `genotype_II_vs_I` will be the comparison between Tg and KO across treatment and placebo,

> `condition_B_vs_A` will be the comparison between treatment vs control for genotype I,

> `genotypeII.conditionB` will be the interaction effect of WT and treatment across genotype III vs genotype I.

> `genotypeIII.conditionB` will be the interaction effect of WT and treatment across genotype III vs genotype I.

I am interested to get the list of genes that are affected by both gene X and treatment ? How can I better model all the effects like treatment and genotype effect on Tg/WT vs KO ?

Answer: So you're interested in condition B vs A, taking care of interaction between genotype II vs I (KO over WT) ? Then it should be `genotypeII.conditionB` result.
Question: I'm trying to convert hundreds of variant positions [found here][1] to vcf for downstream analyses and cannot find a good way to do this. I found [jannovar][2] but that only take variants from .c (coding) and .n (non-coding) positions whereas portal.gdc.cancer.gov produces .g (genome) positions.

An example in hgvs format would be

17:g.7674180C>A
17:g.7675997G>T
17:g.7676257G>A
17:g.7676088G>C
17:g.7676215G>A
17:g.7676152delC
17:g.7676381C>A
17:g.7670712delG
17:g.7670716C>G
17:g.7676264_7676265insA

Any help would be much appreciated.

[1]: https://portal.gdc.cancer.gov/exploration?facetTab=mutations&filters=%7B%22op%22%3A%22and%22%2C%22content%22%3A%5B%7B%22op%22%3A%22in%22%2C%22content%22%3A%7B%22field%22%3A%22ssms.consequence.transcript.annotation.vep_impact%22%2C%22value%22%3A%5B%22HIGH%22%5D%7D%7D%2C%7B%22op%22%3A%22in%22%2C%22content%22%3A%7B%22field%22%3A%22ssms.ssm_id%22%2C%22value%22%3A%5B%22set_id%3AAWoQcB6Ycg-xDiQUTMVg%22%5D%7D%7D%5D%7D&searchTableTab=mutations&ssmsTable_offset=600&ssmsTable_size=100
[2]: https://doc-openbio.readthedocs.io/projects/jannovar/en/master/index.html

Answer: Well, after a bunch of wasted time VEP was the winner. All I had to do was paste my list of hgvs variants to http://uswest.ensembl.org/Homo_sapiens/Tools/VEP and there is a download as VCF option.
Question: Hi,

This might be another naive question. I have RNA-Seq data from 2 different types of samples (pain vs no-pain), and I want to **identify novel lncRNA/eRNA in pain**. I also have DHS/ATAC data for that gives me accessible regions in the DNA, some of which would be regulatory.

One way to do this is to use Trinity (using all samples) for de-novo transcript assembly, then use Salmon **on this assembled transcriptome** to find transcripts that are differentially expressed, and then among those, the transcripts that don't have prior annotation would be 'novel' transcripts. Cufflinks can also be used for this.

Another way would be to just directly run Salmon on the samples, using the **hg38 transcriptome**, and find differentially expressed transcripts, some of which could be novel.

I want to know if doing this via the first method will lead to *more* novel transcripts, and are there any merits to the first method over the second, in general.

Thank you!

Answer: I would not use Trinity on human samples. Trinity is mainly usable when you do **not** have a reference genome.

Since you have a reference genome **and** a well annotated transcriptome you are much better of using tools which perform guided de-novo transcript reconstruction. The most well known tools for that are Cufflinks and StringTie. I have written more about the detailed steps and considerations [here][1].

[1]: https://bioconductor.org/packages/devel/bioc/vignettes/IsoformSwitchAnalyzeR/inst/doc/IsoformSwitchAnalyzeR.html#what-quantification-tools-should-i-use
Question: I have a BAM file that has 33X whole genome coverage.

I would like to create a new BAM file from the original that has 10X genome coverage.

Is this the correct method?

samtools view -h -s 0.3 in.bam >10x.bam

Also what methods do you use for fast calculation of whole genome coverage? GATK is bloody slow.

Answer: Yes, the `-s 0.3` will result in an average 10x coverage. One simple method to get total coverage is to simply parse the output of `samtools depth`.
Question: Hi all!!
I have a large vcf file and I want to create a subset one according to #CHROM field with a txt file (a list that contains #CHROM IDs of interest). I would like to keep the headers and the vcf format.
Any ideas of how to do that?
Thanks a lot! :)

Answer: https://www.biostars.org/p/201603/
https://www.biostars.org/p/46331/
https://bioinformatics.stackexchange.com/questions/3401/how-to-subset-a-vcf-by-chromosome-and-keep-the-header
Question: I have a dataframe (**trip**) that contains a column (**SNP**). It looks like this (but longer, and it has 192 levels):

SNP
C[T->C]T
C[G->C]A
G[A->C]A
C[T->C]C
C[C->A]G
T[G->A]C
...

I want to pattern match and replace on the following criteria:

gsub("G->T", "C->A", trip)
gsub("G->C", "C->G", trip)
gsub("G->A", "C->T", trip)
gsub("A->T", "T->A", trip)
gsub("A->G", "T->C", trip)
gsub("A->C", "T->G", trip)

but ALSO, if one of the patterns listed above is found, I want the string in which it's contained have additional substitutions applied. Namely:

if ((grep(G->T|G->C|G->C|A->T|A->G|A->C), trip$SNP)==TRUE){
substr(trip$SNP, 1,1) <- tr /ATCG/TAGC/; #incompatible perl syntax?
substr(trip$SNP, 8,8) <- tr /ATCG/TAGC/;
}

As in, if any of these patterns--G->T, G->C, G->C, A->T, A->G, or A->C--is found in a string in trip$SNP, replace the 1st and 8th characters in that string according to this regex: tr /ATCG/TAGC/;

Desired output, with changes highlighted:

SNP
C[T->C]T
C[G->C]A
G[A->C]A
C[T->C]C
C[C->A]G
T[G->A]C

to:

SNP
C[T->C]T
G[C->G]T #<-- changed
C[T->G]T #<-- changed
C[T->C]C
C[C->A]G
A[C->T]G #<-- changed

Is there a more elegant way to do this?

Answer: ## assuming df is your data frame
SNP <- as.character(df$SNP)

SNP
[1] "C[T->C]T" "C[G->C]A" "G[A->C]A" "C[T->C]C" "C[C->A]G" "T[G->A]C"

i <- grep("(A|G)->", SNP)
SNP[i] <- chartr("ACGT", "TGCA", SNP[i])

SNP
[1] "C[T->C]T" "G[C->G]T" "C[T->G]T" "C[T->C]C" "C[C->A]G" "A[C->T]G"
Question: Dear all,

I use a tool called ExomeDepth. 2 months ago I was able to run it and it worked well. Today I decided to analyse different samples and try to count reads:

library(ExomeDepth)
data(exons.hg19)
bams <- c(long vector of bams)
my.counts <- getBamCounts(bed.frame = exons.hg19, bam.files = bams, include.chr = F, referenceFasta = "/path/to/fasta/hg19.fa")

It worked before and it is exactly what is written in the manual: https://cran.r-project.org/web/packages/ExomeDepth/vignettes/ExomeDepth-vignette.pdf . I have R-3.6.1 and ExomeDepth 1.12 (freshly installed).

The error:

Error: RangedData objects are deprecated and the rownames setter for
RangedData objects is now defunct. Please migrate your code to use
GRanges or GRangesList objects instead. See IMPORTANT NOTE in
?RangedData

(obviously the important note says that it is deprecated)

Is there a way to fix it manually or I should wait for the update of the package?

(This is not a crucial thing, I run this tool for benchmarking, so if I won't be able to run it I'll just write it like it is, but if smb has thoughts how to fix it - I'd do this)

Answer: You could try installing from vplagnol github repository (using DevTools) to see if the dependency has been updated (GRanges)

[https://github.com/vplagnol/ExomeDepth][1]

install.packages("devtools")

library(devtools)
install_github("vplagnol/ExomeDepth")

If it hasn't then you could try raising an issue on the ExomeDepth repository.

[1]: https://github.com/vplagnol/ExomeDepth
Question: Dear Friends,

I am working with histone proteins in my new project. I am trying to determine the genes interacting with HDAC11 gene. ultimately looking for genes various pathways leads to cancer. is there any good webservers or databases I can look for sorting out the genes for this study?

Thanks for reply

Answer: Hi sana,

Please have a look in " string" database. Here you need to chose specific organism to see the interacting gene.

http://string-db.org/