cannin/biostars_qa huggingface dataset:
Question: HI,

I had a single fasta file with 60,000 bacterial genomes, I'd like to extract the entire fasta sequence which has a keyword **Psychro** in its header ( Ex - *Psychrobacter* - 18982 genomes).

I'm aware we can extract a subseq using the seqtk

seqtk subseq test.fa test.txt

But I want to extract entire fasta sequence by providing the fasta headers in the test.txt

I provided example fasta header names (~18982 totally genomes I want to extract)

>NZ_CAJHBU010000049.1 **Psychrobacter vallis** isolate Psychrobacter vallis CMS39, whole genome shotgun sequence
>NZ_CAJHBM010000029.1 **Psychrobacter sp. JCM** 18903 isolate Psychrobacter sp. JCM18903, whole genome shotgun sequence
>NZ_CAJHBB010000047.1 **Psychrobacter sanguinis** isolate Psychrobacter sanguinis 13983, whole genome shotgun sequence

>NC_007204.1 **Psychrobacter arcticus** 273-4, complete sequence
>NC_007969.1 **Psychrobacter cryohalolentis K5**, complete sequence
>NC_007968.1 **Psychrobacter cryohalolentis K5** plasmid 1, complete sequence
>NC_008709.1 **Psychromonas ingrahamii 37**, complete sequence
>NC_020802.1 **Psychromonas sp. CNPT3,** complete sequence
>NC_018721.1 **Psychroflexus torquis ATCC** 700755, complete sequence

Suggestions please!

Answer: If you have the headers as shown above you can use the answers here: https://www.biostars.org/p/329774/#329776
Question: Hi,

I am conducting a survival analysis using the gene expression values (converted to z-scores) from a certain gene and event free survival time for 249 patients. I am a newbie to this, so please excuse my lack of expertise (this is also my second question on here, so sorry for any formatting errors, I'm still figuring this out).

Here is my code:

> finaldata$SurvObj <- with(finaldata, Surv(efst, status == 2))

> finaldata$SurvObj

> res.cox1 <- coxph(SurvObj ~ Zscore, data = finaldata)

> summary(res.cox1)

And here is the output:

Call:
coxph(formula = SurvObj ~ Zscore, data = finaldata)
n= 247, number of events= 140
(2 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
Zscore 0.5878 1.8001 0.2118 2.776 0.0055
---
Signif. codes: 0 0.001 0.01 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
Zscore 1.8 0.5555 1.189 2.726
Concordance= 0.561 (se = 0.025 )
Likelihood ratio test= 7.92 on 1 df, p=0.005
Wald test = 7.71 on 1 df, p=0.006
Score (logrank) test = 7.74 on 1 df, p=0.005

efst is the event free survival time for patients, the status is whether they are dead or alive (1 is alive, 2 is dead), and Zscore is the Z-Score value for the expression of the gene in question for each of the patients (to understand whether the expression is relatively high or low).

However, I find it highly improbable that I am receiving such p-values and that leads me to think that I have definitely done something incorrectly since this is my first time trying this out.

Please let me know whether you are able to find any errors I have made/how to fix them, or any help at all with how to better interpret this data. I appreciate it!

Answer: Nothing seems out of the ordinary to me. The upper- and lower- confidence intervals of the hazard ratio are even both above 1.0; however, the difference between them is somewhat large, which may indicate heteroskedasticity, outliers, or some other minor problem. With no other information about your data source and processing steps, there's not much else that we can say.

Thanks.

Kevin
Question: Hi all, another noob question from me.

I finally finished that FASTA Splitter from my last post, which I will update soon.
My new task is to Chunk a Fasta File into well... chunks still in Python too.

FileInput = input('What is full address of the file you need to work with? ' )
FileOutPut = input('What is the address of the Directory you want saving to?')
OrganismOI = input('What organism are you looking at? ')
ChunkSize = int(input("What is the chunk size (in kb) you want? "))

Chunk = ChunkSize*1000
Counter = 0
FileName = FileOutPut + '/' + OrganismOI + str(Counter) + '.fa''

FILE = open(FileInput, 'r')
File = FILE.read()
FileSize = len(File)
print(str(FileSize) + ' is the length of the file in bp.')
#This also includes Fasta TitleNumber
print('Final Number of Files to be expected is ' + str(FileSize / Chunk) + '.')

for chunk in File:
Save = open(FileName, 'w')
Save.write(File[:Chunk])
Save.close()

For the life of me, as soon as it comes to loops, I die inside. No matter how much I read or practise it just isn't going in.
This script is supposed to go through the original FASTA file, and every chunk (ChunkSize, being user defined) spit out a file containing that chunk and called OrganismOI + counter (Supposed to count file numbers). Once again and always, help is very much appreciated.
If I can then add a co-ordinate system ill be set for not.
If anyone has any books or materials that could feasibly help me out learning this stuff then please send it my why. Would the BioStar handbook be a sensible purchase?

Again thanks everyone.

Answer: Hi,
try this [link](https://stackoverflow.com/questions/312443/how-do-you-split-a-list-into-evenly-sized-chunks). (just load the genome as string and pass with to one of the implementations of `chunker` from the link).

If you have **very** low RAM or **very** long genome then you could do something like:

f = open(bare_fasta, 'r'):
l = f.readline()
to_write = []
c = 0
i = 0
threshold = 50000
while l:
to_write.append(l)
c += len(l)-1
if c > threshold:
with open('fout{}'.format(i), 'w') as o:
o.write(''.join(to_write))
i += 1
to_write = []
c = 0
l = f.readline()
with open('fout{}'.format(i), 'w') as o:
o.write(''.join(to_write))

This will split your original fast into chunks by whole lines just after the total number of nucleotides exceeded defined threshold. Be aware that this expects FASTA file with many lines. NOT one long line.
Question: Hi all -

Does anyone have ideas for how to detect whether a particular genome has undergone a whole genome duplication? I have an assembly file (scaffolds) from a de-novo genome assembly. I have identified some duplicated proteins, but after reassembling with an assembler for highly heterozygous genomes, most of these duplicated genes disappear. I'm not sure if this actually answers the question regarding whole genome duplication though. Any thoughts are appreciated!

Thank you.

Answer: One common and relatively simple method of looking for evidence of a whole genome duplication is to estimate the rate of synonymous substitutions (Ks) between duplicated gene copies and then plot the distribution of Ks and look for peaks in the distribution. The logic of this is that if many genes are duplicated at once (which obviously occurs during a WGD), there should be a relatively large number of paralogous gene copies with similar Ks values, but if genes are simply being duplicated at some low underlying rate, you would not expect peaks in the distribution.

There are many papers that discuss this type of method (and its limitations). Here is recent one:

https://academic.oup.com/gbe/article/10/11/2882/5100828
Question: Hi,

This is a pretty general question I know, but I'm wondering why I might be seeing a statistically significant peak from a GWAS in one population and nothing in another. Statistical power to detect an association is roughly the same in both populations.

Hope somebody can give me some insight.

Sally.

Answer: The locus is real, but only has an effect in one population. You could hypothesize that this is due to interaction with other alleles associated with that population, or due to an environmental difference between the populations, or due to differences in the condition. This has been seen several times (though no one has adequate explanation that I'm aware of) in common diseases such as various cancers. Alternately, the effect is real in both populations but the effect is so modest that other environmental or genetic factors hide the signal in one population. That's a different situation from true interaction.

The locus, despite apparent statistical significance, is in fact only appearing by chance. That's why you need replication for these studies.

There is a systematic difference in genotyping technique between populations, e.g. different platforms, different DNA extraction methods, hybridization at different genome cores, etc. This difference is affecting this locus, either creating a new signal or removing one selectively.

The significant probe is flawed somehow: it cross-hybridizes, it has weak signal that is noisy, etc. This gives a false positive result in one population.

Have a look at the raw genotype data to see if the calls are credible in both populations. Look for biases associated with population. Look for systematic clinical/macroscopic differences in your phenotype between populations For example, if this were breast cancer, is one population African-derived and another European-derived? If so, not only do you have different allele frequencies as Jorge points out, you also have different predispositions to ER-positive vs ER-negative cancer (probably due in part to the genetic architecture, but also due to environmental differences). If both populations were Caucasian, is there an ascertainment bias due to different medical systems diagnosing patients differently? Etc...
Question: Hello, I have bulk rna-seq data that consist of two conditions (control, treated) tested across three time points (day 0, 1 and 3). We already observed by simple pair-wise (intra-condition) EdgeR that a big number of genes exhibit a change in expression between each two **adjacent** time points (day 1 vs 0 , then day 3 vs 1). We want to find significant differences INTER-conditions i.e genes whose change is significantly different in Treated vs Control at day 1 vs 0 , then day 3 vs 1. How could I do this with EdgeR? I have run EdgeR before but I do not know how to test both influence of treatment and time in one single design formula, and make sure that adjacent time points be tested (and NOT 3 vs 0 for exemple). I have 3 biological replicates by time-point and condition. Thanks for any help!

Answer: You could use a full-factorial design by combining the time variable and treatment variable like (simplified exampe):

sample	treatment	time	factorial

sample1	case	t0	case_t0
sample2	control	t0	control_t0
sample3	case	t1	case_t1
sample4	control	t1	control_t1
sample5	case	t3	case_t3
sample6	control	t3	control_t3

and then design as `~factorial` followed by making the desired contrasts. Contrasts could be something like:

`(case_t1 - control_t1) - (case_t0 - control_t0)`

That would test for genes where the difference between case and control is differential between time points. For example a modest change between case and control at t0 that is much greater in t1 would could come out as significant.
This is called an interaction. See e.g. section 3.3 in the manual for more details and examples, https://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf
Question: Hi everyone,

I want to count the reads per allele using pysam.The script I have pasted here is only for a variant,but I am doing for many variants in the main script.I want to filter out the reads that overlap,so that I am not double counting them.Basically I want to check the difference in the count of reads per allele with vs without overlap.This can not be done using pysam functionality,so I am trying to do something without it.Below is my script:

#! /usr/bin/python
import sys
import random
import pysam
samfile = pysam.Samfile("f1.bam", "rb")
out = open('outfile', 'w')
f= open("f2.vcf", "r")
ref_counter=0
alt_counter=0
ref='G'
alt='A'
total=0
for pileupcolumn in samfile.pileup( 'chr1',881626,881628):
if int(pileupcolumn.pos) == 881626:
if int(pileupcolumn.n) <= 300:
for pileupread in pileupcolumn.pileups:
if str(pileupread.alignment.seq[pileupread.qpos]) == str(ref):
ref_counter = float(ref_counter) + 1
elif str(pileupread.alignment.seq[pileupread.qpos]) == str(alt):
alt_counter = float(alt_counter) + 1
total=ref_counter+alt_counter
print(str(ref) + '\t' + str(ref_counter) +'\t' + str(alt) + '\t' + str(alt_counter) + '\t' + str(total) + '\n')

Any suggestions?

Answer: I assume you have paired-end reads and don't want to double count things if both read of a pair cover some position of interest. The simplest method would be to just add read names to a list as you process them. You can then check that before processing and skip an alignment if the read name is already in there.

for pileupcolumn in samfile.pileup('chr1',881626,881628):
found = []
...
for pileupread in pileupcolumn.pileups:
try:
found.index(pileupread.alignment.qname)
except ValueError :
if str(pileupread.alignment.seq[pileupread.qpos]) == str(ref):
...
found.append(pileupread.alignment.qname)
...

That's the general idea at least.
Question: I have tried to draw a Venn diagram like below, but I faced with Error :

> 'ERROR [2019-04-01 17:10:26] Impossible: a2 <- area2 - a9 - a10 - a19
> - a20 - a21 - a11 - a28 - a29 - a31 - a22 - a30 - a26 - a25 - a24 - a14 produces negative area. Error in draw.quintuple.venn(area1 =
> 327498, area2 = 30526, area3 = 639, : Impossible: a2 <- area2 - a9 -
> a10 - a19 - a20 - a21 - a11 - a28 - a29 - a31 - a22 - a30 - a26 - a25
> - a24 - a14 produces negative area'.

grid.newpage()
venn.plot <- draw.quintuple.venn(
area1 = 327498,
area2 = 30526,
area3 = 639,
area4 = 2000,
area5 = 49925,
n12 = 34944,
n13 = 130,
n14 = 305,
n15 = 31838,
n23 = 350,
n24 = 1014,
n25 = 19339,
n34 = 64,
n35 = 214,
n45 = 932,
n123 = 57,
n124 = 267,
n125 = 12847,
n134 = 3,
n135 = 31,
n145 = 179,
n234 = 41,
n235 = 186,
n245 = 706,
n345 = 30,
n1234 = 3,
n1235 = 22,
n1245 = 189,
n1345 = 3,
n2345 = 27,
n12345 = 3,
category = c("A", "B", "C", "D", "F"),
fill = c("dodgerblue", "goldenrod1", "darkorange1", "seagreen3", "orchid3"),
cat.col = c("dodgerblue", "goldenrod1", "darkorange1", "seagreen3", "orchid3"),
cat.cex = 2,
margin = 0.05,
cex = c(1.5, 1.5, 1.5, 1.5, 1.5, 1, 0.8, 1, 0.8, 1, 0.8, 1, 0.8, 1, 0.8,
1, 0.55, 1, 0.55, 1, 0.55, 1, 0.55, 1, 0.55, 1, 1, 1, 1, 1, 1.5),
ind = TRUE
)

Thanks for your help in advance!

Answer: why not simply use [Vennerable][1] with method "AWFE"

to be clear, the error indicates a miscalculation in the intersect areas. If you want, you can use the function `Venn` from the mentioned library and use these numbers. I would simply plot the resulted venn using `plot(..., type = "AWFE")`.

hth

[1]: https://r-forge.r-project.org/scm/viewvc.php/*checkout*/pkg/Vennerable/inst/doc/Venn.pdf?revision=58&root=vennerable
Question: I have list of regions ( chr:start-end) for which I would like to build two consensus haplotypes. I could not use the FastaAlternateReferenceMaker from GATK and it ignores the phasing information with in the region of interest. My question is more detailed in this image:

https://www.dropbox.com/s/hkg4eff9ba8zg37/Screen%20Shot%202015-12-26%20at%2013.57.49.png?dl=0

I could try GATK HaplotypeCaller but I have already called SNPs using freebayes and would like to know if there is any existing method/tool to do this. There are tools that constructs haplotypes but they look complicated and are for different purpose. I am looking for something simple. Otherwise I will end up writing a dirty script, which may not work in all possible scenarios.

Answer: just a quick suggestion: use ReadBackedPhasing ( https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_tools_walkers_phasing_ReadBackedPhasing.php ) to generate a VCF with the phasing information. Use this phasing to produce two VCFs for each strand, and then go back to FastaAlternateReferenceMaker ?
Question: Hi,
I'm trying to learn pybedtools using the intersect and am very new to it.

I want to use the -wb option but can't quite work out how to do this with pybedtools.

At the moment I have:

import pybedtools

a = pybedtools.BedTool(en.bed)
b = pybedtools.BedTool('annot.bed')

a.intersect(b)

If I do:

a.intersect(b, -wb )

I get the error

> 'NameError: name 'wb' is not defined'

Answer: I got an answer from the developer [https://github.com/daler/pybedtools/issues/266][1]

[1]: https://github.com/daler/pybedtools/issues/266
Question: I am looking to detect if there is GSTM1 deletion in CEU population from 1000 genome project. Can anyone advice where to start? I have tried to extract CNVs from 1000 genome vcf file using vcftools but it outputted only 5 deletions none of which were on chromosome 1. I was thinking about downloading sequence files but I am not sure what to download? Can anyone advice?

Answer: On the VCF file for chr1 from 1000 genomes, since you already know the gene as well as the population, why not just do the following:

1. GATK's `SelectVariants` to pick only the genotypes corresponding to your population (you can get the genotype to population info from the 1000 genomes website - I can help you get there if you wanna do this)
2. vcftools with `--chr`, `--start` and `--end` to pick the gene region

That should give you a manageable region that you could use UNIX utils with, even eyeball.

Or, you could zoom in to the region on the UCSC genome browser, enable 1000 genomes track in full mode and gain a visual perspective on it.

Also, 1000 genomes has VCF files per chromosome. Maybe you ran the CNV extractor on the wrong file?
Question: Hi guys,

does anyone know is there a way to access ATCC database to get all the info about the cell lines, like tissue of origin, species etc? Or is there another db that will cover all that? I tried hyperCLDB but still a lot of cell lines are not covered in there.

And I know that I can search the ATCC website but that works with 2 or 3 cell lines and not with a bunch of them.

Will be grateful for any tips!
Kasia

Answer: I don't know how to delete this post and I found what I was looking for, so I will share:

https://web.expasy.org/cellosaurus/
This database seems to be the most comprehensive one from the ones I looked at and it comes in different file formats (xml, txt or obo).

Also, all the information about the cell lines used by Encode can be found here: https://genome.ucsc.edu/ENCODE/cellTypes.html
Question: Simple question.... Has anyone produced a plugin for adding [ExAC][1] allele frequencies to Variant Effect Predictor (VEP) annotation?

[1]: http://exac.broadinstitute.org/

Answer: As the data are available as a VCF file, you can just use the custom annotations feature of the VEP to pull in the sites and any frequency data you want:

```
curl -O ftp://ftp.broadinstitute.org/pub/ExAC_release/release0.3/ExAC.r0.3.sites.vep.vcf.gz
curl -O ftp://ftp.broadinstitute.org/pub/ExAC_release/release0.3/ExAC.r0.3.sites.vep.vcf.gz.tbi
perl variant_effect_predictor.pl -i example_GRCh37.vcf -cache -custom ExAC.r0.3.sites.vep.vcf.gz,ExAC,vcf,exact,0,AF
```

Here I've added the AF field (it will appear as ExAC_AF in the VEP output), but you could add any of the fields in the INFO block of the VCF by listing them at the end of the `--custom` string. See http://www.ensembl.org/info/docs/tools/vep/script/vep_custom.html for more info.

You could get a list of fields by doing (this gets the VCF header)

tabix -h ExAC.r0.3.sites.vep.vcf.gz 1:1-1

There's a caveat here in that the VCF carries the per-population allele counts but not the frequencies, so to derive the frequencies you'd need to do a little math(s). It's probably worth writing a plugin to deal with these (and a very simple one it would be too, see the dbNSFP plugin http://www.ensembl.info/ecode/dbnsfp/ for an example of one that uses tabix to retrieve data from a file), but in the short term the above solution should at least get you part of the way.
Question: Hello,

When using NCBI's interface I can exclude species, genus, even phyla. Is there a way to do that in command line blast?

I need to exclude one species from my search. I see this can somehow be done with a GI list, but the example they give is not very understandable.

Adrian

Answer: Use NCBI's taxonomy browser and link from there to genbank GI list of all nucleotides or proteins. Download that and use it in the command line blast, for example from here to here, then "send to" -> "file" -> "GI list"
Question: Hi,

I have a table and a vector, e.g.,

Table =

A B C
s1 0 1 2
s2 1 2 1
s3 2 0 1

Vector:

vec <- c(0.01, 0.087, 0.34)

I wish to iterate over the table to isolate pairs of rows and then use the vector to fit a 1st order linear regression model `(Y~A*B)`. My problem is that I wish to initiate a matrix outside the main for loop and append the variable 'p_value' to that matrix after each iteration.

The code to isolate the p-values is given below:

for (rows in 1:nrow(Table)) {
li_row <- Table[rows,] # Extract a single row
row <- unlist(li_row) # Convert to vector

for (other in 1:nrow(Table[2, ])) {
li_row_2 <- Table[other,] # Extract a different single row
row_2 <- unlist(li_row_2)
lin_reg <- lm(vec ~ row*row2) # Call 1st order linear model
summ <- summary(lin_reg) # Store results summary
f_stat <- summ$fstatistic # Extract F-statistic data
p_value <- fstat[1] # Extract p-value [**The value to append to matrix**]
}
}
}}

In short, upon each iteration, I wish to add the 'p-value' to a matrix. I know a priori how big the matrix will be and also know what labels to attach to rows and columns.

Thanks,
S ;-)

Answer: Quick and dirty answer; concatenate all your p-values into a single vector as you progress

vec <- c(vec, p_value)

and at the end, convert to a matrix

matrix(c(0, 1, 2, 1, 2, 1, 2, 0, 1), nrow = 3)

Better answer; don't iteratively append, but rather map a function down the rows of your table and use R's built-in vectorisation to create a matrix.

A toy example below applies an anonymous function to each row `y` of matrix `x`, returning a new matrix. The function adds the row `y` to another row, also from `x`, but chosen randomly (without replacement).

x <- matrix(c(1, 2, 3, 10, 10, 10, 4, 5, 6), nrow = 3)
apply(x, 2, function(y) y + x[sample(1:nrow(x), 1), TRUE])
Question: Hi,

I have a matrix of protein expression that has some missing values. I know its MNAR (missing not at random). However, I am having trouble using the impute.QRILC command to impute the matrix as I get the same matrix as a result. Would someone be able to me understand why the matrix is not imputing? My code is the following:

y_qrilc <- impute.QRILC(norm_prot)[[1]]

`norm_prot` is my protein expression matrix

Any help is appreciated.

Answer: Multiple imputation is recommended for MNAR data (see the link to the paper below). You can use the mice package in r to do this.

[Multiple Imputation Article][1]

[1]: https://bmcmedresmethodol.biomedcentral.com/articles/10.1186/s12874-017-0442-1
Question: I'm using Bowtie2 for mapping reads to human genome. I'm using Bowtie2 first time, so what are the best parameters to set if we are using these mapping results to further use to detect SNPs? bowtie2 has also given some preset options. which is more appropriate for finding SNPs from alignment results.

Which one is more appropriate local or end-to-end mode? if my purpose of alignment is SNP detection?

Answer: Just use one of the presets, likely either --very-sensitive or --very-sensitive-local. If you've not trimmed the reads (or just trimmed off adapters), then you should use local alignment. If, however, you've adapter & quality trimmed the reads, then using end-to-end alignment would seem more relevant. Realistically speaking, you should get very similar results (from variant calling) regardless. In general, local alignment will sometimes give a bit better results, since you can still get alignments in cases of structural rearrangements or variants at the end/start of reads.
Question: Hi all,

I have a vcf file looks like

chr22 31600460 rs2413035 T C,TCTTCATCAGCCTTCTCCTCCTCTTCCTC .

I'm trying to expand the two variants into two seperate rows like

chr22 31600460 rs2413035 T C .
chr22 31600460 rs2413035 T TCTTCATCAGCCTTCTCCTCCTCTTCCTC .

Of course I want other fields(INFO, FORMAT) to be expanded accordingly. "genotype" has to be expanded as well.

Does anyone have any idea or comments?

Thanks.

Answer: You can use bcftools norm for splitting multi-allelic variants, try something like this:

```
bcftools norm -Ov -m-any input.vcf > output.vcf
```