cannin/biostars_qa huggingface dataset:
Question: I am trying to download a number of GSE datasets with GEOquery using codes that looks like the following

library(GEOquery)
library(Biobase)
library(limma)

gset <- getGEO("GSE54884", GSEMatrix =TRUE, getGPL=FALSE)

and I am receiving this error:

Error in open.connection(x, "rb") : HTTP error 403.

Here is my session information:

> sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows >= 8 (build 9200)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] parallel stats graphics grDevices utils datasets methods base

other attached packages:
[1] GEOquery_2.50.5 Biobase_2.42.0 BiocGenerics_0.28.0

This seems to work for other GSE datasets such as:

gse <- getGEO("GSE781",GSEMatrix=FALSE)

My understanding is that HTTP error 403 is an access issue, but I thought all GEO Datasets are open access. Any ideas about what I might be doing wrong?

Answer: It looks like the problem has been resolved.
Question: Dear lazyweb,

my question is somehow related to some previous questions asked here: ( http://biostar.stackexchange.com/questions/788 , http://biostar.stackexchange.com/questions/31, http://biostar.stackexchange.com/questions/613 , etc...).

I'm looking for a table that would summarize the pro's and the con's of each tool for NGS (mappers and/or snp callers), the algorithm used, if the quality of the fastq, is taken into account, how the gaps & the indels are handled, a publication, etc...
... e.g.: something that would look like the following (almost empty) page on wikipedia: http://en.wikipedia.org/wiki/List_of_sequence_alignment_software#Short-Read_Sequence_Alignment

Thanks,

Pierre

UPDATE: @Julien_f suggested this link: http://lh3lh3.users.sourceforge.net/NGSalign.shtml

Update2: another cool resource cited in Marc's paper: http://seqanswers.com/wiki/Software/list

Answer: Although there are over 30 aligners, I guess 90% of researchers only choose from the list below:

Illumina: Bowtie, BWA, Mosaik and Novoalign (Stampy is also a good one)

SOLiD: Bfast, Bowtie and BWA (perhaps also NovoalignCS?)

454: SSAHA2, LAST, BLAT and Mosaik (I would recommend BWA-SW actually)

Carefully studying their features would not take much time. In addition, to "feel" a tool, we'd better use it rather than just look at the feature table.

As to SNP callers, I guess 90% who do not write in-house tools are using GATK and SAMtools/MAQ. SNVMix is also a good one, especially for cancer sequencing. As there are only a few of them, it is easy to try each.

In addition, usual SNP calling on deep resequencing data is not hard. Many groups use their in-house tools and being highly tuned for specific data, in-house tools may have better performance.

At last, as I am promoting everywhere, apply BAQ implemented in samtools. Doing this will improve the accuracy for all SNP callers so far.
Question: Is there way to combine these samtools commands:

1. `samtools view -h -F 0x10 `

2. `samtools calmd `

I tried this:

samtools view -h -F 0x10 | samtools calmd | samtools view -bS - >

But it doesn't work.

Answer: In bash you need to use `-` to specify that the input is coming from a pipe. You should also add the `-S` parameter to specify that the input is SAM (or follow the suggestion of @swbarnes2 and output BAM from `view` instead):

samtools view -h -F 0x10 | samtools calmd -S - | samtools view -bS - >
Question: Hello every one, can any one suggested a tool or method to get rsID for a number variants I have in vcf file, do I have to manipulate the header for my file? Many thanks for your help in advance

The vcf file I have in this shape:

```
##fileformat=VCFv4.1
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##LeftAlignVariants="analysis_type=LeftAlignVariants bypassFlowAlign=true kmer_len=19 min_var_count=5 short_suffix_match=5 min_indel_size=4 max_hp_length=8 min_var_freq=0.15 min_var_score=10.0 relative_strand_bias=0.8 output_mnv=0 sse_hp_size=0 sse_report_file= target_size=1.0 pref_kmer_max=3 pref_kmer_min=0 pref_delta_max=2 pref_delta_min=0 suff_kmer_max=3 suff_kmer_min=0 suff_delta_max=2 suff_delta_min=0 motif_min_ppv=0.2 generate_flow_position=0 input_file=[] read_buffer_size=null phone_home=STANDARD gatk_key=null read_filter=[] intervals=null excludeIntervals=null interval_set_rule=UNION interval_merging=ALL reference_sequence=/results/referenceLibrary/tmap-f3/hg19/hg19.fasta rodBind=[] nonDeterministicRandomSeed=false downsampling_type=BY_SAMPLE downsample_to_fraction=null downsample_to_coverage=1000 baq=OFF baqGapOpenPenalty=40.0 performanceLog=null useOriginalQualities=false BQSR=null defaultBaseQualities=-1 validation_strictness=SILENT unsafe=null num_threads=1 combined_sample_name= num_cpu_threads=null num_io_threads=null num_bam_file_handles=null read_group_black_list=null pedigree=[] pedigreeString=[] pedigreeValidationType=STRICT allow_intervals_with_unindexed_bam=false logging_level=INFO log_to_file=null help=false variant=(RodBinding name=variant source=/results/analysis/output/Home/MGHBED_602/plugin_out/variantCaller_out/IonXpress_001/small_variants.sorted.vcf) out=org.broadinstitute.sting.gatk.io.stubs.VCFWriterStub NO_HEADER=org.broadinstitute.sting.gatk.io.stubs.VCFWriterStub sites_only=org.broadinstitute.sting.gatk.io.stubs.VCFWriterStub filter_mismatching_base_and_quals=false"
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##fileDate=20140616
##phasing=none
##reference=/results/referenceLibrary/tmap-f3/hg19/hg19.fasta
##reference=file:///results/referenceLibrary/tmap-f3/hg19/hg19.fasta
##source=Torrent Unified Variant Caller (Extension of freeBayes)
#CHROM POS ID REF ALT QUAL FILTER
chr1 65886142 . C G 97.35 PASS
```

Answer: Use [GATK VariantAnnotator][1] with `--dbsnp`

[1]: https://www.broadinstitute.org/gatk/guide/tooldocs/org_broadinstitute_gatk_tools_walkers_annotator_VariantAnnotator.php#VariantAnnotations
Question: I have a PLINK formatted database (bed/bim/fam files) and corresponding
recoded files (hh/ped/map). I am looking for an effective way to extract
samples from this database with specific genotypes. I have looked through the PLINK manual
and found that I can extract set of samples using "--keep" parameter and
extract set of genotypes using "--extract", am wondering if this can be
done in a single step using another parameter or tool.

My input is a list of rsIDs and genotypes; I need to get sample ids and
genotype as output.
INPUT

rs1800562 AA

OUTPUT

Sample1 AA
Sample5 AA
Sample22 AA
...

Is there any option in PLINK to do this or I should use unix 'grep' and/or a custom script to extract data. Suggestions on other computational genomics tools to do similar task is also welcome.

Answer: You can combine both --keep and --extract in a single step, but you're wanting to condition your --keep based on the genotypes you get from your --extract, which PLINK can't do to my knowledge. If you want a single ped file for each snp you could do something like

awk '{print $1}' INPUT > mysnps
plink --bfile data --extract mysnps --tfile mysnps
(some code here to loop through each line of mysnps.tped and pulling out column index when your genotype matches, and write out a list of samples for each snp)
(some code here to run plink --keep for each list of samples)

... but you probably already knew this, and just need an implementation. Sorry this wasn't much help.
Question: Hi! I am actually looking at TCGA level3 RNASeqV2 data. My goal is to look at the DEGs (tumor vs. normal) and I'm looking at LUAD now.

I am using edgeR at the moment since the original rsem paper mentioned that those rsem can be processed by edgeR/ DESeq.

I have a couple questions that I was wondering if anyone might have any suggestion -

1. Does it make sense to include all the tumor samples available, including those that don't have the matching normal samples from the same participant, and analyze for the DEGs? What kind of normalization method would be recommended if I do so? Or can I just use the default normalization of edgeR?

2. I started out looking at only the matched TN and NT samples. Using the above I'm getting 5639 DEGs out of 20531 genes (FDR <=0.05, FC >=2) which seems like a lot? ( even a lot more if I don't use any FC filter)

3. There seems to be various discussion regarding what tools to use:
- http://seqanswers.com/forums/showthread.php?t=28515
- https://groups.google.com/forum/#!topic/rsem-users/H1cswrvvmPs

I wonder if anyone has more experienced in analyzing TCGA dataset has any thought as to whether it is OK to use EdgeR, or should I use some other tools like EBSeq for RNASeqV2 data?

Any suggestion is greatly appreciated. Thanks a lot in advance!

Answer: 1. My personal opinion is that including all the samples would not be good Idea, you would have lot of variance and so would loose many DEGs. Once your matrix is ready go for quantile normalization then you should go for hierarchical clustering of samples with spearman correlation. Choose biggest (should have much more samples than other clusters) clusters (one from each patient and normal). Make these two groups and do DE analysis.
2. Yes, you are getting lot of DEGs , If you want to narrow down, then I have seen some publications where people take FDR threshold 0.01 also. But its okay 1/4 of total genes.
Question: Hello happy bioinformaticians :)

It can be a very simple question but I want to ask that how can I count line (row) of header of VCF ?

I can be done manually but I want to get accurate result.

Thanks,
BG

Answer: Simply count the number of rows starting with a `#`:

grep -c '^#'

If you want to count the number of sequence headers:

grep -c '##contig'

And even more.. if you want to count the number of non-headers:

grep -c '^[^#]'

~~On the command line of course ;)~~
Question: Hi, I have two ChIP-seq bed files of same transcription factor. I want to compare one specific binding motif between these two ChIP-seq binding peaks, i.e. whether the motif is more enriched in the binding sites of one ChIP-seq vs. the other (under different culture condition)?

The plots I want to make is by X axis: "Distance from centre of ChIP region (bp)" and the Y axis: " No. motif / no. binding sites"

Thank you for your answers!

Xiaoyong

Answer: You can use MEME suite for enrichment of a specific motif over the given background (http://meme.nbcr.net/meme/tools/ame).. To plot the enrichment with distance, you will have to bin your input sequences into categories (1kb, 2kb, 3kb from TSS etc.), compute enrichment over control and then plot it in R with bins on X and enrichment score on Y axis..

Hope this helps.
Question: Hey guys,

I was trying to learn genome profiling newly (since a week ago!!!), suddenly few days ago I knew that I can do some of my analysis in galaxy free from my OS and all.. does this mean that I can do my supervisor's tips **all** in galaxy and without typing code in cmd?

Answer: I'll summarize the linked thread by saying "not necessarily." Here are some questions you'll want to ask yourself (or your supervisor). If the answer is "yes," then it works in Galaxy's favor:

- Is the analysis workflow rigid and constant?
- Are the workflow results few in number and small in size?
- Can we avoid making complex or nuanced decisions about the next step at every point in the analysis workflow?
- Is there a Galaxy instance out there (there are others besides the main public one at usegalaxy.org) which has all of the tools that we need?

If Galaxy still seems appropriate after asking these questions, then it would likely be useful to peruse the public workflows in the "Workflows" tab to see if anyone's replicated all or part of what you want to do.

Lastly, if you have specific questions about Galaxy, [it has its own Biostar instance where a lot of really helpful and knowledgeable people hang out][1].

[1]: https://biostar.usegalaxy.org/
Question: Dear all, I have an account on Illumina basespace and my Run is completed.

Now I want to download the entire run and run bcl2fastq conversion myself. I tried the approaches (my system is mac) mentioned [here][1] but I have not been successful to make the download happen.

Everytime I get this error when running executable bs command: ```ERROR: *** stat /Users/.basespace: no such file or directory ***```

Could you please give some input about how I can do it? Thanks

[1]: https://developer.basespace.illumina.com/docs/content/documentation/cli/cli-overview

Answer: Hi,
You must have downloaded the cmdline utility bs.
The first step is to authenticate your executable.
You probably have done this (?)

$ ./bs auth

When you do that, a unique URL is printed on the cmdline, which you then need to copy into your browser. You are then presented with the sign-in portal of Illumina basespace. Provided that you have had created a log-in already, you need to sign-in. Once you do that, a welcome msg gets printed on the cmdline, saying like "Welcome username".

You are then ready to go. If authentication happened successfully, then under your home directory a file would have been created. It would look like -

$ cat /home/user/.basespace/default.cfg
apiServer = https://api.basespace.illumina.com
accessToken = xxxxxxx54343b3xxxxxa5exxx

Then, you can list out the project thats available in your Ilmn BaseSpace dir. for downloading by doing -

$ ./bs list project
+----------------------------------------------------+-----------+-------------+
| Name | Id | TotalSize |
+----------------------------------------------------+-----------+-------------+
| 22.10.04_my_prj_xxxxx1_2022-10-04T09_02_03_1xxx137 | 36xxxxxx4 | 36618128088 |
+----------------------------------------------------+-----------+-------------+

The number under the Id header could then be used to download all fastq under that project. Like -

./bs \
download \
project -i 36xxxxxx4
Question: UPDATE: Partly solved. It works without sorted option and now I do not need this option. However, the original question remains the same - how to calculate coverage from .bam file using `-sorted` option. I added dummy region from MT to the .bed file - it did not work.

Have a problem with coverageBed. Run command:

covBed -d -sorted -b Sample.realigned.dm.recalibrated.bam -a srt.bed

Get output:

```
ERROR: Database file Sample.realigned.dm.recalibrated.bam contains chromosome , but the query file does not.
Please re-reun with the -g option for a genome file.
See documentation for details.
```

How to solve it? (https://www.biostars.org/p/159721/ - was not really helpful). It works without `-sorted` option.

My **bed file** chromosomes:

```
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
X
Y
```

**Bam file**:

```
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
X
Y
MT
```

Answer: Have you actually tried to pass a genome file to coverageBed? As I remember, coverageBed will complain even if chrom names are the same and in the same order but one chrom is missing from one of the two inputs, as for example chrom MT in your case. You can obtain the genome file from UCSC with e.g:

mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -e \
"select chrom, size from hg19.chromInfo" > hg19.genome

Or from the bam file header with:

samtools view -H aln.bam \
| grep -P "@SQ\tSN:" \
| sed 's/@SQ\tSN://' \
| sed 's/\tLN:/\t/' > genome.txt

Then run `coverageBed -g genome.txt ...`
Question: GenCode has a wonderful breakdown of the number of coding/non-coding genes and transcripts here...

http://www.gencodegenes.org/stats.html

Does a similar breakdown exist for RefSeq anywhere? Thanks!

Answer: If you take a look at this [hg19/GRCh37 GFF file from NCBI][1] here is what is in there:

495245 CDS
1 D_loop
33314 Genomic
1 RNase_MRP_RNA
1 RNase_P_RNA
2 SRP_RNA
4 Y_RNA
22 antisense_RNA
19086 cDNA_match
618452 exon
30661 gene
6405 lnc_RNA
49861 mRNA
22621 match
3097 miRNA
31 ncRNA
2046 primary_transcript
23 rRNA
297 region
1 sequence_feature
72 snRNA
393 snoRNA
698 tRNA
1 telomerase_RNA
5732 transcript
3 vault_RNA

For GRCh38 [a similar file can be found here][2] with following statistics:

1413848 CDS
1 D_loop
33893 Genomic
26 Genomic%2CXM/XP/XR
1 RNase_MRP_RNA
1 RNase_P_RNA
2 SRP_RNA
15 V_gene_segment
4 Y_RNA
22 antisense_RNA
13572 cDNA_match
24 centromere
1856502 exon
43504 gene
28117 lnc_RNA
114575 mRNA
22834 match
3038 miRNA
31 ncRNA
2025 primary_transcript
23 rRNA
558 region
304 sequence_feature
62 snRNA
389 snoRNA
629 tRNA
1 telomerase_RNA
16011 transcript
3 vault_RNA

[1]: ftp://ftp.ncbi.nih.gov/genomes/Homo_sapiens/GRCh37.p13_interim_annotation/interim_GRCh37.p13_top_level_2017-01-13.gff3.gz
[2]: ftp://ftp.ncbi.nih.gov/genomes/Homo_sapiens/GRCh38.p10_interim_annotation/interim_GRCh38.p10_top_level_2017-01-13.gff3.gz
Question: Hello,

A colleague pointed out this problem with GATK, which generates indexes in memory on the fly if it doesn't find those on disk:

INFO 14:05:14,023 RMDTrackBuilder - Creating Tribble index in memory for file dbsnp_132.vcf
WARN 14:09:31,798 FSLockWithShared - WARNING: Unable to lock file dbsnp_132.vcf.idx (could not open read/write file channel)
WARN 14:09:31,798 RMDTrackBuilder - Unable to write to dbsnp_132.vcf.idx for the index file, creating index in memory only

I've read and gone through these two resources:

http://www.broadinstitute.org/gsa/wiki/index.php/Tribble

http://code.google.com/p/tribble/

But I found it a bit awkward having to write my own java class just to generate an index file from a .vcf :-!

VCFTools does not seem to have that functionality either, at least according to their documentation:

http://vcftools.sourceforge.net/docs.html

Anybody knows how to generate those .vcf.idx files in a more straightforward way ?

Thanks in advance !

PD: I just couldn't refrain myself from including this link too :)

http://en.wikipedia.org/wiki/The_Trouble_With_Tribbles

Answer: Hi
I found IGVTools, best for VCF indexing.

igvtools can be run from the command line or IGV itself (File>Run igvtools...) After launching, choose the Index command and browse to your .vcf file. The index file (.idx) will be created in the same directory as the .vcf file.

Thanks!
Jitendra
Question: I'm using edgeR for DE analysis of RNASeq data. I've been using manual filtering suggested in the EdgeR manual till now.

I followed genefilter package and tried to chose a theta that has the optimal number of rejections. Basically this is what I did: run edgeR once for the data set, then use the p-value results for genefilter calculations and plots.

```
dge <- DGEList(counts=counts[,!is.na(group)], group=group[!is.na(group)])
d <- dge
d <- calcNormFactors(d)
d <- estimateDisp(d, trend="none", robust=TRUE)
et <- exactTest(d, pair = levels(d$samples$group))
dgeTable = topTags(et, n=nrow(d))

res = data.frame(
filterstat = rowSums(d$counts[rownames(dgeTable),]),
pvalue = as.data.frame(dgeTable)$PValue,
row.names = rownames(dgeTable))

rejBH = filtered_R(alpha=alpha, filter=res$filterstat, test=res$pvalue, theta=theta, method="BH")

theta_chosen = theta[max(which(rejBH==max(rejBH)))]
use = res$filterstat > quantile(res$filterstat, probs=theta_chosen)
dgeFilt=dge[use,]
d <- dgeFilt
```

I used rowSums as filtering criteria, and the after I filter out the genes - in this case I got theta = 0.19, I ran the edgeR exact test again:

```
### rerun the model
d <- dgeFilt
d <- calcNormFactors(d)
d <- estimateDisp(d, trend="none", robust=TRUE)

et <- exactTest(d, pair = levels(d$samples$group))
```

My question is - even though the previous calculated max(rejBH) - the max number of rejections is 102 , using the same alpha (0.1) I'm only getting 60 DE genes after applying the filter. Actually, without applying any filter, the number of DE genes is 75 (this matches rejBH at 0%).

- Why is the new DE genes number not 102?
- Why am I getting less DE genes even though from the rejection plots it looks like I'll have the max number of rejections at theta = 0.19?

Am I doing anything wrong here or did I misunderstand something?

Thanks in advance for any suggestion!

Answer: You only need to run exactTest() once and you don't need to subset d. Once you have your filter threshold, you just readjust the pvalues of the included rows in res.

BTW, the reason the numbers differ when you subset d is because you now have a different background distribution used for shrinkage.
Question: Hi,

Using [biomaRt][1] package in R, I have extracted the Transcription Start/Stop Sites of protein coding genes in human genome. So, I have a huge list of TSS. The next step is to get sequences of 1000 bp region upstream and downstream of retrieved TSS coordinates. I figured out that there isn't a function in biomaRt to download sequence upstream and downstream of the TSS.
I do appreciate if you could make a suggestion to me.

[1]: https://www.bioconductor.org/packages/devel/bioc/html/biomaRt.html

Answer: You can use the appropiate BSGenome from Bioconductor (probably [BSgenome.Hsapiens.UCSC.hg19](http://bioconductor.org/packages/BSgenome.Hsapiens.UCSC.hg19/) or [BSgenome.Hsapiens.UCSC.hg38](http://bioconductor.org/packages/BSgenome.Hsapiens.UCSC.hg38/). You should:

1. Transform your list of list of TSS to a GRanges object

2. Get the upstrem or downstream region with [flank](https://bioconductor.org/packages/release/bioc/vignettes/GenomicRanges/inst/doc/GenomicRangesIntroduction.html)

3. Use [getSeq](https://bioconductor.org/packages/devel/bioc/manuals/BSgenome/man/BSgenome.pdf) to retrieve the actual sequence from BSgenome

With this approach you'll only download the sequence once, when downloading the BSgenome package instead retrieveing it from the Biomart server every time you rerun the script.
Question: Is there any software available for simulating a bacterial population which contains similar strains?

Answer: I am currently doing some simulations where I need bacterial populations at varying levels of identity to each other, so I wrote a tool to do that in the BBMap package. For example, using E.coli:

mutate.sh in=e_coli.fasta out=mutant.fasta id=99
randomreads.sh ref=mutant.fasta out=reads.fq.gz reads=5m length=150 paired adderrors

That will create a mutant version of E.coli with 99% identity to the original, and then generate 5 million simulated read pairs from the new genome. You can repeat this multiple times; each mutant will be different.