cannin/biostars_qa huggingface dataset:
Question: Hello All,

Let's say that some whole genome sample was sequenced with a coverage of 30x.
As far as i'm aware, this means that, with respect to the reference genomes' nucleotides,
the data represents each nucleotide 30 times on average.

Let's also say that the tissue sample was heterozygous for some loci,
where the frequency of the two alleles are both 0.5.
Does this mean that coverage for each of these locations are, in effect
15x?
I.e if you aligned the data (and it aligned correctly), you would
expect to see ~15 reads with allele 1 and ~15 with allele 2.

N.B

I ask because I am trying to make a simulated cancer genomics dataset.
For this I am using ART, and have "mutated" the hg19.fa file, by introducing
some point mutations.
This mutated file with represent one haploid set, whilst the non-mutated hg19.fa
file will represent the other haploid set; this should add realistic point mutations, which
are usually heterozygous in nature.

I then plan to sequence at 30x, so, I was going to run ART for each file at 15x and
then combine to get 30x. Any thoughts?

Thanks,
Izaak

Answer: Theoretically the coverage per allele would indeed be total coverage/2 (for diploid genomes). However, that's more often not exactly 50%.
Question: Hi,

I wonder why for calculating gene lengths, they calculate the median of the transcripts in goseq manual? Its under section 5.3. I suppose it should be sum rather than median. The manual can be found here.

http://www.bioconductor.org/packages/release/bioc/vignettes/goseq/inst/doc/goseq.pdf

Cheers!

Answer: The goal is to determine if there's a bias by gene length. In order to do this, one needs to derive some sort of gene length measure. There are a couple ways to do that:

Union gene model: The total non-redundant exonic length of a gene

Estimated length: Derived from using expectation maximization, where you then have an estimate of the expected gene length within each sample

Median transcript length: What's used here, which is just the median of the annotated transcript lengths.

If you actually summed the transcripts, as you suggest, then you'd get odd results for genes with many isoforms. The result also wouldn't match any biologically plausible length (e.g., if a gene has 20 isoforms, each ~1kb then your method would yield a length of 20kb rather than a more reasonable ~1kb estimate).
Question: Hi,

I am trying to match the location of the list of string like this

library(seqinr)

at <- ("ATATATAT")
s1 <-ifelse(at[8]=="T"||"A" && at[7]=="A"||"T" &&
at[6]=="T"||"A",5,
ifelse(at[2]=="T"||"A" && at[4]=="A"||"T" &&
at[1]=="T"||"A",'1','0'
))
s1

It works fine only for one sequence. I tried it in a for loop but getting error like

invalid 'x' type in 'x && y'

Any help is much appreciated
Thanks

Answer: To simplify your example, condition: if any 2nd or 4th position in every sequence has A or T, then TRUE.

# example data
x <- c("AAGTA",
"AAGTA",
"AAGTA",
"ACGAA")

# in this example all TRUE
all(substr(x, 2, 2) %in% c("A", "T") | substr(x, 4, 4) %in% c("A", "T"))
# [1] TRUE

If this is not the solution you are looking for, then please provide example input and expected output, clearly.
Question: I am trying to post query to a webserver : http://www.imtech.res.in/raghava/antibp/submit.html

but I am getting an error

Traceback (most recent call last):
File "crawler.py", line 4, in
conn = httplib.HTTPConnection("http://www.imtech.res.in/raghava/antibp/submit.html")
File "/usr/lib/python2.7/httplib.py", line 704, in __init__
self._set_hostport(host, port)
File "/usr/lib/python2.7/httplib.py", line 732, in _set_hostport
raise InvalidURL("nonnumeric port: '%s'" % host[i+1:])
httplib.InvalidURL: nonnumeric port: '//www.imtech.res.in/raghava/antibp/submit.html'

The python script is shown below:

import httplib, urllib
params = urllib.urlencode({'seqname':'GICACRRRFCPNSERFSGYCRVNGARYVRCCSRR','format':'Amino acid sequence in single letter code', 'terminus':'N-terminus', 'method':'svm', 'svm_th':'0', 'type': 'Submit'})
headers = {"Content-type": "application/x-www-form-urlencoded", "Accept": "text/plain"}
conn = httplib.HTTPConnection("http://www.imtech.res.in/raghava/antibp/submit.html")
conn.request("POST", "", params, headers)
response = conn.getresponse()
print response.status, response.reason
data = response.read()
conn.close()

What could be the problem? Thank you.

Answer: That's how I would do it, with the disclaimer that I'm no expert in querying web pages and I don't know anything about the server in question:

python
import mechanize

br = mechanize.Browser()
br.set_handle_robots(False)
br.open("http://www.imtech.res.in/raghava/antibp/submit.html")
br.select_form(nr = 0)

## See what is available on this web page:
for f in br.forms():
print f

#
#
# )>
#
#
#
#
#
#
# =Submit) (readonly)>
# =)>>

## Input your sequence and parameters:
br['seqname']= 'myseq'
br['seq']= 'GICACRRRFCPNSERFSGYCRVNGARYVRCCSRR'
br['format']= ['nformat']
br['terminus']= ['1']
br['svm_th']= '0'

## Sumbit and collect results:
res= br.submit()
html= res.read()

Now `html` is string in html format that you could parse with an html parser or something else. The relevant bit in `html` should look like:

Antibacterial Activiy

GICACRRRFCPNSER	1	1.975	YES

GYCRVNGARYVRCCS	18	1.051	YES

ICACRRRFCPNSERF	2	1.001	YES

...
Question: Hi, all

I use awk to extract rows from a text file:

awk 'NR==FNR{vals[$1];next} ($1) in vals' indiv.txt file1.txt > new_file1.txt

But how can I use the same code for multiple files ?

I can use :

awk 'NR==FNR{vals[$1];next} ($1) in vals' indiv.txt file1.txt > new_file1.txt
awk 'NR==FNR{vals[$1];next} ($1) in vals' indiv.txt file2.txt > new_file2.txt
awk 'NR==FNR{vals[$1];next} ($1) in vals' indiv.txt file3.txt > new_file3.txt
awk 'NR==FNR{vals[$1];next} ($1) in vals' indiv.txt file4.txt > new_file4.txt
awk 'NR==FNR{vals[$1];next} ($1) in vals' indiv.txt file5.txt > new_file5.txt
. . . . . . . . . . . . . . . .
awk 'NR==FNR{vals[$1];next} ($1) in vals' indiv.txt file29.txt > new_file29.txt

But there is some ways to do it automaticaly ? for example a loop or file*.txt > new_file*.txt.

Thanks!

Abdel

Answer: There is a variable called `FILENAME` which provide the current filename `awk` reads in. You can take this to create a filename `awk` should print its output.

$ awk 'NR==FNR{vals[$1];next} ($1) in vals { print $0 >> "new_"FILENAME}' indiv.txt file1.txt file2.txt file3.txt

fin swimmer
Question: Why is it so difficult to make things in ggplot2 , i like the way it helps in customisation but the curve is steep nevertheless

Here is my sample dataframe

df <- gene HSC CMP
ENSG00000158292.6 1.8102636 2.456869
ENSG00000162496.6 2.6796705 6.203838
ENSG00000117115.10 3.4509115 5.555739
ENSG00000159423.14 3.6809277 5.063446
ENSG00000053372.4 5.7089974 6.851090

If i have plot a boxplot i can simply write this `boxplot(df[,-1],col=c("red","blue"))`

I get a boxplot but when im trying with ggplot2 im having difficult time

ex <- melt(df, id.vars=c("HSC", "CMP"))
ggplot(data = ex,
aes(x = CMP, y = HSC)) +
geom_boxplot()

I get a single boxplot what i want is i get a box plot for HSC and CMP which i got when i use simple base R boxplot .

Any help or suggestion would be highly appreciated with my ggplot2 code

Answer: ex = melt(df, id.vars="gene")
ggplot(ex, aes(x=variable, y=value)) + geom_boxplot()

Your `melt()` command produced nonsensical output.
Question: I can easily view the **chromHMM tracks** for **Roadmap project** but not able to download the bed files from the UCSC genome browser. I want the **bed files** for **H3K27ac marks** for a couple of tissues like **Heart Left Ventricle**, **Fetal Heart** etc. I am attaching a screenshot of the browser with Heart Left Ventricle chromHMM track. Also, [this is the link][1] for the same at UCSC Genome Browser just in case the screenshot does not upload. The track is second from top and names chromHMM tracks from Roadmap.

[![< image not found >][3]][2]

Does anyone know how to download these tracks as bed files?

[1]: https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&position=chr10%3A180000-187400&hgsid=392137799_1PzPhmDS4ZXa40AFMunFuaQTgkjR
[2]: https://www.dropbox.com/s/bmjeuhp200vkklc/chromHMM_Roadmap_H3k27ac_HeartLeftVentricle.png?dl=0
[3]: https://www.dropbox.com/s/bmjeuhp200vkklc/chromHMM_Roadmap_H3k27ac_HeartLeftVentricle.png?dl=0

Answer: I had emailed Ting Wang (of the **Roadmap Epigenomics Project**) and he promptly replied back addressing my question. Here is the answer to this question, in his own words:

> All the original chromHMM tracks are here for you to download: https://sites.google.com/site/anshulkundaje/projects/epigenomeroadmap
>
> All Roadmap/ENCODE data are available through the WashU EpiGenome Browser (http://epigenomegateway.wustl.edu/). If you right-click on any track and select "information", you will usually see a url for data source.
Question: Hi everyone,

I'm working on Copy Number data from TCGA. I download "Gene Level Copy Number Variation" using TCGABiolinks R package and the following code:

library(TCGAbiolinks)

query_cnv <- GDCquery(project = "TCGA-KICH",
data.category = "Copy Number Variation",
data.type = "Gene Level Copy Number Scores")
GDCdownload(query_cnv)
data <- GDCprepare(query_cnv)

Everything works great. I get a nice dataframe with first three columns being: "Gene.Symbol" / "Gene.ID" / "Cytoband".
To facilitate the analysis and being able to merge data from other sources such as RNASeq, I tried to convert Ensembl gene ids contained in Gene.Symbol in Hugo Symbol using BioMart.

library(biomaRt)
mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
genes <- gsub(".\\.","\\1",data$Gene.Symbol)
geneIDs <- getBM(filters = "ensembl_gene_id", attributes = c("ensembl_gene_id","hgnc_symbol"), values = genes, mart = mart)

However, over 19729 different Ensembl ID, I only get 3269 match.

What is surprizing is that according GDC docs (https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/CNV_Pipeline/), this dataset should contain CNV associated to each gene, so I would expect a little bit more match according coding regions.

When I tried to search the description of Ensembl ID not found by Biomart. I get zero answers from both Ensembl and NCBI. (Example: "ENSG000000081221" "ENSG000000081314" "ENSG000000676014" "ENSG000000783616" "ENSG000000788015"); So, it's like these ID did not exist in any database. So, where they are coming from ?

Did I miss something ? Is it normal to have few genes coding for proteins in this kind of datasets ? Should I process differently for the analysis of such data ?

Any suggestions or comments will be really helpful.

Answer: The IDs as you present them do not exist.

The problem is this line of code, which is not doing what [I believe] you believe it's doing:

genes <- gsub(".\\.","\\1",data$Gene.Symbol)

You need to remove the final digits from each Ensembl ID after the dot. This works:

library(TCGAbiolinks)
query_cnv <- GDCquery(project = "TCGA-KICH",
data.category = "Copy Number Variation",
data.type = "Gene Level Copy Number Scores")
GDCdownload(query_cnv)
data <- GDCprepare(query_cnv)
data <- data.frame(data)

ens <- sub('\\.[0-9]*$', '', data$Gene.Symbol)

require(org.Hs.eg.db)
ens_to_symbol <- mapIds(
org.Hs.eg.db,
keys = ens,
column = 'SYMBOL',
keytype = 'ENSEMBL')
head(ens_to_symbol)
ENSG00000008128 ENSG00000008130 ENSG00000067606 ENSG00000078369 ENSG00000078808
"CDK11A" "NADK" "PRKCZ" "GNB1" "SDF4"
ENSG00000107404
"DVL1"

library(biomaRt)
mart <- useDataset('hsapiens_gene_ensembl', useMart('ensembl'))
ens_to_symbol_biomart <- getBM(
filters = 'ensembl_gene_id',
attributes = c('ensembl_gene_id', 'hgnc_symbol'),
values = ens,
mart = mart)
ens_to_symbol_biomart <- merge(
x = as.data.frame(ens),
y = ens_to_symbol_biomart ,
by.y = 'ensembl_gene_id',
all.x = TRUE,
by.x = 'ens')

head(ens_to_symbol_biomart)
ens hgnc_symbol
1 ENSG00000000003 TSPAN6
2 ENSG00000000005 TNMD
3 ENSG00000000419 DPM1
4 ENSG00000000457 SCYL3
5 ENSG00000000460 C1orf112
6 ENSG00000000938 FGR

Kevin

===============
-----------------
Question: Dear Biostars,

For those familiar with IBM LSF, I am wondering if someone knows how to configure a job, via the job itself or a config file, how to handle stdout when a job enters SSUSP states. As far as I understand my issue, a job will restart from the beginning after entering SSUSP but the stdout is appended and not overwritten. Below is an example of a submission. This happens for a variety of functions regardless of output type (compressed or plain text file)

bsub -R 'span[hosts=1]' -N -o stdout.file myFunction ...

Answer: In my specific case, passing `-oo` instead of `-o` was a solution because this overwrites instead of appends the output, however I was not able to address the original issue, will update if I find a work-around...

https://scicomp.ethz.ch/wiki/LSF_mini_reference
Question: Hello everyone,

I want to make a GSEA using the online tool of the Broad Institute.
To do that I generally use FPKM or RPKM matrix and I was wondering if you think it's ok to use the normalized read file that I get as output from DESeq2

thank you

Answer: If you have already run DEseq2 it would be easier for you to use GSEA with the "preranked" mode, giving as input a metric based on the logFC and the pvalues that you got as result of your differential expression.

It is pretty well explained here : http://crazyhottommy.blogspot.com/2016/08/gene-set-enrichment-analysis-gsea.html
Question: Hello everyone,

I would like to download somatic SNP data from the TCGA. But if I have a look at the data matrix right here, there are two color codes "Tumor, matched normal" and "Normal, matched Tumor". I looked up the online guide and the Getting Started with the Data Matrix guide.

They explain it like

TN (Tumor, matched normal) – Data for a tumor tissue for which matched normal tissue exists.

NT (Normal, matched tumor) – Data for normal tissue for which matched tumor tissue exists.

But where is the difference?

May some of you are more experienced using TCGA then me.

With all the best,

Mario

Answer: There is no difference between TN and NT for somatic mutations, because tumors and normals are paired up for somatic variant calling. It only makes sense for data generated separately for tumors/normals... like RNA-seq or methylation assays. The data matrix doesn't have a very intuitive interface because they tried to generalize the filtering UI for different data types.

I'd recommend downloading the somatic mutations (MAF files) directly from the [TCGA DCC][1] or from [Firehose][2]. Read more about TCGA MAF files [here][3].

You can also use the [firehose_get][4] script to download latest available MAFs that are fed into Firehose. Here's how to do that:

```
wget http://gdac.broadinstitute.org/runs/code/firehose_get_latest.zip
unzip firehose_get_latest.zip
./firehose_get -b -only Mutation_Packager_Calls data latest
```

To list the other kinds of data (like `Mutation_Packager_Calls` in the command above) look them up [here][5].

[1]: https://tcga-data.nci.nih.gov/tcgafiles/ftp_auth/distro_ftpusers/anonymous/tumor/
[2]: http://gdac.broadinstitute.org/runs/stddata__latest/samples_report/index.html
[3]: https://www.biostars.org/p/69222/
[4]: https://confluence.broadinstitute.org/display/GDAC/Download
[5]: http://gdac.broadinstitute.org/runs/stddata__latest/data/
Question: Hi Guys,

I am using pileup in Rsamtool from R package. I need to get the total "count" for each position. For example, for chr11, position 1643082, there are a total of 1 count of A, 3 counts of C and 18 counts of T. I want them in this order:

chr11:1643082 A:1, C:3, T:18

and so forth for position chr11, 47663948 to be shown in the result. There are thousands of different positions and I need to calculate the counts for each nucleotide.Thanks!

`table1`

```r
seqnames pos strand nucleotide count which_label
1 chr11 1643082 + A 1 chr11:1643082-1643082
2 chr11 1643082 - C 3 chr11:1643082-1643082
3 chr11 1643082 + T 15 chr11:1643082-1643082
4 chr11 1643082 - T 3 chr11:1643082-1643082
5 chr11 47663948 + C 16 chr11:47663948-47663948
6 chr11 47663948 - C 11 chr11:47663948-47663948
7 chr11 47663948 + T 2 chr11:47663948-47663948
```

Answer: You can do this easily with `data.table` package.

Something like this:

```r
df[, sum(count), by=list(seqnames,pos,nucleotide)]
```
Question: Hi,

I am going to get some data from plasmid sequencing to identify SNPs on the plasmids. What it is done on the lab is the following:

- The plasmids are purified by size
- We amplify the plasmids using the phi29 polimerase. The polimerase will go trough the plasmid multiple times. Hence we get the same sequence concatenated multiple times. **My question is related to this step**
- We sequence it using Oxford Nanopore

My question is:

On step two, I wrote that "we get the same sequence concatenated multiple times". For me, this a potential source of information to correct the base calls prior to aligning them to the reference. Since you have the same sequence multiple times ( concatemers). What I would like to know is if there is a tool that uses the information of the concatemers to improve the base calls. I have try to find some methods but found none. I know that PacBio has a similar flavour using the "circular consensus calling", But I have not found any methodological explanation.

thanks!

Answer: Thanks to this [answer][1] in bioinformatics stackexchange, there is a tool doing exactly what you want. The tool is called [C3P0a][2]

[1]: https://bioinformatics.stackexchange.com/questions/5397/error-correction-within-the-long-read/5406#5406
[2]: https://github.com/rvolden/C3POa
Question: Hello everyone,

I want to download all rbcL sequences for a species list. I don't have the Genbank ID, I just want to use a list of species name, then download all the rbcL genes, which were deposited in Genbank.
I appreciate if anyone can give a suggestion.

Best,

Lingyun

Answer: I resolved the problem. The NCBIminer is pretty well
Question: I have 600,000 + unique protein sequences that I want to blast against themselves. For example, I want to run sequence 1 against sequence 1, sequence 2 against sequence 2, ... sequence 3453 against sequence 3453 etc.

I have local blast installed on my computer and blastdb is prepared for the protein sequences. Problem is, each sequence has to be blast against the whole database when I will only retain the hit from a unique sequence against itself. Needless to say, I need a much more efficient method! Any idea?

Answer: ```
curl -ks "ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz" | gunzip -c |\
awk '/^>/{printf("\n%s\n",$0);next;} {printf("%s",$0);} END{printf("\n");}' | grep -v '^$' |\
while read T
do
echo "$T" > prot.fa
read S
echo "$S" >> prot.fa
blastp -query prot.fa -subject prot.fa
done
```
Question: Hello,

I've watched a video from YouTube about OMIM. Then I get that OMIM already had more than 3000 disease-associated genes till last year. Now I want to get all of those genes from OMIM to do some experiments.

I think maybe I should use "Gene Map" of OMIM, but I still don't know how to search to get all disease-associated genes.

Do you have any good advice or guide?

Thank you in advance,

naulty

Answer: You can download MIM ids and corresponding GeneID from [here][1].

[1]: http://www.omim.org/downloads
Question: I have the vector:

stri = c("Stage II", "Stage IVA", "Stage IIIB", "ACB" )

and I want to remove the A and B in "Stage IVA" ,"Stage IIIB", but not in "ACB",

I want to get

c("Stage II", "Stage IV", "Stage III", "ACB" )

Is there any method by using gsub?

Answer: Fortunately, I find the way to solve it.

Here is the answer:

gsub("(^Stage.*)(A|B)","\\\\1", stri, perl = TRUE)
Question: Hi there!

I have a list of paired fastq. I want to filtrate human read from each of these pairs using BBmap tool.
here is a functionnal snakemake rule i wrote:

#Removing human reads from fastq
rule clean_fastq:
message:
"Removing human reads from fastq."
input:
unzip_fastq_R1 = rules.fastq_unzip.output.unzip_fastq_R1,
unzip_fastq_R2 = rules.fastq_unzip.output.unzip_fastq_R2,
BBMAP = rules.get_BBmap.output.BBMAP ,
index = rules.create_index.output.index
output:
R1_cleaned = result_repository + "FASTQ_CLEANED/{sample}_R1_cleaned.fastq",
R2_cleaned = result_repository + "FASTQ_CLEANED/{sample}_R2_cleaned.fastq"
params:
path_human= result_repository + "FASTQ_CLEANED/"
shell:
"""
{input.BBMAP} in1={input.unzip_fastq_R1} in2={input.unzip_fastq_R2} \
basename={rules.clean_fastq.params.path_human}{wildcards.sample}_%.fastq outu1={output.R1_cleaned} outu2={output.R2_cleaned} \
path=temp/
"""
Each jobs have a cost of about 20Go RAM and the probleme is that I can only have 32 Go available.
I dont know if this is possible for snakemake to execute all jobs from a same rule in a queue, to avoid this memory problem.

If not, I probably should check an other tool to process these fastq. Any ideas? (except bmtagger, I had too much problem with it haha)
What would you suggest?

Thx,

Hadrien

Answer: I think you can use the [resources](https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#resources) directive in combination with the `--resources` command line option. E.g., your rule `clen_fastq` could be:

rule clean_fastq:
resources:
mem_gb= 20,
input:
...

then:

snakemake -j 10 --resources mem_gb=32 ...

This will run at most 10 jobs at a time using at most 32 GB of memory (of "mem_gb" in fact) which means `clean_fastq` will run one a time.
Question: Hello.
I am new to differential gene expression analysis.
I was trying the tutorial of griffith test data provided in BIOSTAR HANDBOOK.
I have done the alignment using HISAT2 and when I am doing quantification using the

> featureCounts

I am getting the GZIP ERROR:-2 as follows:

> GZIP ERROR:-2

the data is paired ended so I am using the same command given in the manual of subhead package.

>featureCounts -p -a annotation.gtf -t exon -g gene id -o counts.txt mapping results PE.bam

Also, the path and working directory all are fine, still, I don't know why this error is coming up.
Please suggest as I am stuck at this step and not able to proceed.

Thanks in advance!

Answer: It seems the error encounter due to `BAM` file input with some bad blocks. Did you check if the validity of the `BAM` file that is generated correctly and not damaged/corrupted? I'd suggest trying a couple of things:

1. Could you try re-running the command with `SAM` file as input in `featureCounts`?

2. Could you try running the command on a different BAM file (maybe from publically available data?)

Have a look at a discussion forum regarding similar error: https://groups.google.com/forum/#!topic/subread/S4smWRfBNPM

Hope this helps!
Question: I have a multiple sequence alignment represented as a "mulit-fasta." I want to remove any column that contains a gap.

Thought I would ask before writing my own.

Answer: Trimal is very versatile. Try with "-gt 0.0"
http://trimal.cgenomics.org