cannin/biostars_qa huggingface dataset:
Question: Hello Biostars Community,

General question(s) here, specifically in regards to sequencing platforms and some related questions on sequencing, as well?

1. Could using a newer/the newest computational genome/annotation (for example, presently, Ensembl 107 or the newest Gencode version) adversely effect the actual truth of what was sequenced?

2. When sequencing is done through an Illumina machine or other big name company machines, are those sequencing platforms *completely independent* from the genome or DNA/cDNA being sequencing?

3. What happens if, for example, "famous gene ABC" and "low-profile gene XYZ" are found to have different 3' and/or 5' ends by some new discovery, in Illumina, would adapters still link to them to perform those bridge PCR reactions on the flow cell lanes, or would it be that all the data published before on "famous gene ABC" and "low-profile gene XYZ" should be revisited? Or is it like question #2, "completely independent" - are even adapters independent from genes?

I was doing some reading, and I guess the gene sequence really only matters for probe-based sequencing (chips and arrays?). Hopefully this question could be a good resource for others?

Thank you in advance.
- Pratik

Answer: Excellent responses by GenoMax to wonderful questions by Pratik. I echo what GenoMax says.

IMHO, the platforms ar einherent to capture technologies as well. NOT necessarily or everytime one would intend to use the same samples with different sequencing platforms/technologies unless the lab is rich. Even otherwise, I wouldn't recommend that. Taken together, the need of the our is to emply machine learning heuristics and predict such outcomes on a large scale. Could completely independent be dependent can give some interesting reasons for us to explore.

Just my two cents
Prash
Question: What are currently the best tools for lifting over genome coordinates from one custom genome build to another custom genome build.

I found this 4 year old post that mentions **UCSC liftOver** as the top tool but that seems to be limited to model organisms as you need the chain file from UCSC . I can't find anywhere how to create these chain files your self:
https://www.biostars.org/p/65558/

Is it possible to create these UCSC liftOver chain files by your self for custom genome builds?

Are there other tools with the same quality and functionality as UCSC liftOver that support creating these "chain" files yourself for custom reference genomes? How difficult and time / resource consuming is it to build these chain files?

Do these chain files and the lift over process take strand into account? If I have a VCF with an A/T SNP positions do they flip it to T/A SNP position if the new reference genome has a genome region flipped to the different DNA strand?

Does lift over support INDEL variants?

Answer: As mentioned elsewhere, our https://github.com/wurmlab/flo tool can do a decent chunk of the job
Question: When I try this [vcftools_results]$ /home/X/vcftools_0.1.6/perl/vcf-validator sample1.flt.vcf.gz

it gives an error.

Can't locate Vcf.pm in @INC (@INC contains: /share/lib64/perl5/5.8.8 /share/lib64/perl5 /share/lib/perl5/5.10.0/x86_64-linux-thread-multi-ld /share/lib/perl5/5.10.0

/share/lib/perl5/5.8.8 /share/lib/perl5 /usr/lib/perl5/5.8.8 /usr/lib/perl5

/share/lib/perl5/5.10.0/x86_64-linux-thread-multi-ld /share/lib/perl5/5.10.0

/share/lib/perl5/site_perl/5.10.0/x86_64-linux-thread-multi-ld

/share/lib/perl5/site_perl/5.10.0 /share/lib/perl5/site_perl/5.8.8

/share/lib/perl5/site_perl .) at /home/X/vcftools_0.1.6/perl/vcf-validator line 9.

BEGIN failed--compilation aborted at /home/X/vcftools_0.1.6/perl/vcf-validator line 9.

how can I fix it?

Answer: Hi Deniz

Have you set the environment variable PERL5LIB as it is defined in VCF Tools installation ?

It is as simple as calling:

export PERL5LIB=/path/to/your/installation/perl

Regards.
Question: Hi,

I am attempting to make some heat maps with some RNA-Seq data. I used Salmon for my quantification and followed [this][1] link to import my data into R so I can use Deseq2. So I am just going through the Deseq2 little tutorial [here][2]. I am trying to follow their instruction on how to create some heat maps. I believe my issue is with the row names. Someone asked a similar question to mine [here][3] but the issue was never resolved. So in that post, it was suggest that the row names were set to the default 1,2,3,4,5,6,7,8, which is true. This is the result I get when I run `rownames(df)`

[1] "1" "2" "3" "4" "5" "6" "7" "8"

So I trued to change the column names by running the following

rownames(df) <- c("adp1","adp2","adp3","adp4","fed1","fed2","fed3","fed4")
dds@colData@rownames <- c("adp1","adp2","adp3","adp4","fed1","fed2","fed3","fed4")

However, I am still getting the same error

select <- order(rowMeans(counts(dds,normalized=TRUE)),
decreasing=TRUE)[1:20]
df <- as.data.frame(colData(dds)[,c("condition")])
pheatmap(assay(ntd)[select,], cluster_rows=FALSE, show_rownames=FALSE,
cluster_cols=FALSE, annotation_col=df)

`Error in check.length("fill") : 'gpar' element 'fill' must not be length 0`I

think my attempts to change the row names were point less? I think the issues goes way back to important my Salmon data into R? I could be mistaken but when I went through their import tutorial, the samples were labeled as sample1, sample2, sample3, etc.

Any help would be amazing!

[1]: https://bioconductor.org/packages/devel/bioc/vignettes/tximport/inst/doc/tximport.html#salmon-sailfish
[2]: http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#heatmap-of-the-count-matrix
[3]: https://www.biostars.org/p/281586/

Answer: You can just do this:

cdata <- colData(dds)

pheatmap(assay(ntd),
cluster_rows = FALSE,
show_rownames = FALSE,
cluster_cols = FALSE,
annotation_col = as.data.frame(cdata))

...or:

pheatmap(assay(ntd),
cluster_rows = FALSE,
show_rownames = FALSE,
cluster_cols = FALSE,
annotation_col = as.data.frame(cdata[,"condition"], row.names=rownames(cdata)))

You do not even have to create the *df* object.

Kevin

***Edit (October 24, 2018:**
the key is that the rownames of the object passed to `annotation_col` have to match the colnames of the object being clustered (here, `assay(ntd)`)*
Question: My lab has started using the Haloplex technology to capture the regions of interest.

With this technology most of the reads have the same start/end (I've found that removing the duplicates would reduce the coverage to ~1 base). The reads are grouped in a vertical 'cluster' that can contains more than 1000 reads.

Have you ever used this technology ? I'm currently using samtools mpileup with the option
"-A" ( count anomalous read pairs) and "-d 8000" ( max per-BAM depth).

I'm afraid there is a kind of bias with haloplex: is it safe to use samtools or Gatk to call the variants ?

Pierre

Answer: In general, any capture-based technology will have bias. In your case, you're worried that the the bias will be amplified in the detection step.

I would use as many methods as possible to call these alleles, and see what the differences are. With coverage this high, there should be little variation in results.

Furthermore, simulation could provide an avenue to assessing potential bias in the case of haploplex resequencing. In simulated tests, I have observed that samtools and GATK have lower sensitivity than our caller (freebayes) at high depths, but all perform roughly the same at lower coverage. It's not entirely clear to me why this is, but I suspect default parameter selection may play a role. I haven't evaluated other callers, and it's possible that others may perform even better in high-depth contexts.
Question: I am looking for any two strains of S. pneumoniae for which I can find their raw output sequencer data in FASTQ format and for which I can find their reference genomes (since I plan to align the raw data in the FASTQ format to the reference genomes using Bowtie 2). I know that I could make artificial FASTQ files using a generator but I need to work with data that came out of a real sequencer.

I tried looking in NCBI, Ensembl, and ArrayExpress and I did find reference genomes for some strains such as [R6][1] and [others][2]. However, finding the raw data straight out of a sequencer (FASTQ format) seems a hard task.

In summary, where can I find any two strains of S. pneumoniae for which can I find the aforementioned raw data (I assume that there must be popular strains for which this is more easily found than for others) along with their reference genomes?

[1]: http://www.streppneumoniae.com/getntseq.asp
[2]: http://strepneumo-sybil.igs.umaryland.edu/cgi-bin/current/shared/index.cgi?site=strepneumo

Answer: You can get FASTQ and FASTA [here][1].

Most strains will provide the raw data associated with the study. For example: http://www.ebi.ac.uk/ena/data/view/ERP000241

[1]: http://www.ebi.ac.uk/ena/data/search?query=streptococcus%20pneumoniae
Question: I should start by saying I have solved this problem, but I feel like my code is ugly and overkill.

I am trying to sort a dataframe on a column containing human chromosomes:

chr1,chr2...chrY, chrX.

the problem is:

chr1, is followed by chr10.

What tricks do you use to deal with this problem IN the [R] environment?

Answer: Hi Zev,

what you are looking for is: mixedsort {gtools}

Order or Sort strings with embedded numbers so that the numbers are in the correct order

package 'gtools' was built under R version 2.13.2

n<- c('chr1','chr21','chr13','chr4','chr10')
> n
[1] "chr1" "chr21" "chr13" "chr4" "chr10"

> mixedsort(n)
[1] "chr1" "chr4" "chr10" "chr13" "chr21"

I hope this helps.
Question: I am trying to characterize the depth of coverage for SNPs/INDELs/other variants in my BAM/SAM file. Using samtools mpileup generates a nice VCF that gives me information on this, but it's not quite complete. Entering the "-t AD" flag gives me the depth of the alternate alleles compared to the reference. That is what I want, though the actual depth is found in the column after the flag is indicated and it is colon delimited. It is readable, but not terribly concise.

However, I am trying to check specific positions along the genome for specific alternate alleles. Is there a convenient way to parse the VCF for the AD flag values by position and allele without complicated and format-specific if/else statements?

I can use bedtools intersect or a dozen other methods to reduce the file to specific positions, already, but I'm hoping there's a method that doesn't require piping a handful of tools together.

EDIT:

It would be straightforward to write an if statement of the kind:

if ALT column is T then return AD and if AD is > 20 return "likely to be target organism strain"
However, over entire genomes and many different variant calls, this won't be so trivial. I have a file that contains a list of positions along the genome and alleles that I'm specifically looking for. I guess I should better word my question to ask for some software that performs that level of analysis. I need to compare allele identity at a given position to a file containing many positions and alleles of interest and then return the coverage of those alleles in the sample that match the ones in the comparison file. After that, I need to review the depth to determine the significance of the variant - which is arbitrary.

Answer: > It would be straightforward to write an if statement of the kind
> if ALT column is T then return AD and if AD is > 20 return "likely to be target organism strain"

using vcffilterjs https://github.com/lindenb/jvarkit/wiki/VCFFilterJS (will add NOT_TARGET_ORGANISM_STRAIN in the column FILTER)

java -jar dist/vcffilterjs.jar -e 'function accept(v) {var foundT=false;var alts= v.getAlternateAlleles(); for(i=0;i< alts.size();++i) {if(alts.get(i).getDisplayString().equals("T")) {foundT=true;break;}} if(!foundT) return false; var g0=v.getGenotype(0); if(!g0.hasAD()) return false; var AD=g0.getAD(); for(var i in AD) if(AD[i]>20) return true; return false;} accept(variant);' -F NOT_TARGET_ORGANISM_STRAIN input.vcf
Question: I have downloaded GnomAD vcf file and want to calculate the minor allele frequency of each variant in that file. But I am not quite sure how to do that. I have an INFO column which has AC,AF, AN in there. Does AF refers to minor allele frequency for that variant? If not then how can I calculate the minor allele frequency. I've given below three variants from the file just to show you what a gnomAD file looks like

`chr1	905315	.	C	A	293.0	PASS	AC=1;AN=143276;AF=6.97954e-06;variant_type=snv...	chr1	905311	906011	chr1:905311-906011	63	.	905431	905432	0,0,0	2	70,528	0,172`

`chr1	905315	rs1489213119	C	T	1432.0	PASS	AC=1;AN=143276;AF=6.97954e-06;variant_type=snv...	chr1	905311	906011	chr1:905311-906011	63	.	905431	905432	0,0,0	2	70,528	0,172`

I am sorry if this is a basic question. I am new to bioinformatics and haven't been able to find a solution yet.

Answer: Yes AF is the MAF here.

It bears mentioning that sometimes the variant might not necessarily be the less common allele in all populations, and the major allele might even confer something weakly deleterious. But in these two cases, they only saw one individual with this mutation.
Question: I am using tabix to retrieve information in a specific genomic region. need not the hole information but only let say second and the third column.

My file is constructed as follows:

```
##INFO=
##INFO=
#SNPID CHROM POS GENE A1 A2 STUDIES CASES CONTROLS
```

So, ideally I would like to have only SNPID, CHROM, POS, CASES.

Is there any possibility to retrieve only this information?

Answer: tabix file.tsv.gz chr12:345-678 | grep -v "##" | cut -f 1,2,3,7
Question: I ran into an unexpected feature of blastn.
After extracting some gene sequences from a genome, creating a blast database and blasting back to the reference, lot of extracted genes are not found in the blast result, while they are certainly there in the genome (as they were extracted from there)
What can be the cause?
Some details (I cannot upload the whole files):

My command is:

blastn -query GCF_000005845.2_ASM584v2_genomic.fna -db MG1655_genes -outfmt 6

The fasta file is downloaded from https://www.ncbi.nlm.nih.gov/genome/167?genome_assembly_id=161521.

The database is gained from extracting the feature table belonging to the assembly above.

A missing gene from the blast is e.g. `aaaD`. However, if blasting only this gene, it is naturally found.

Answer: After doing some research, I found the answer for my question. The value for `max_target_seqs` is 500 by default.
If raising `max_target_seqs` to some irrationally high value, all genes are shown.

Hence I used

blastn -query GCF_000005845.2_ASM584v2_genomic.fna -db MG1655_genes -outfmt 6 -max_target_seqs 100000000

to obtain all genes.
Question: Hi,

I analyse micro array datasets to get gene expression data.

I want to filter NCBI datasets by sample type.

For example I want to get only tumor/cancer sample in [GSE19826][1].

I look at all the samples one by one and select them manually.

Is there any simple way or tool to get it or I have to look manually?

![enter image description here][2]

[1]: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE19826
[2]: https://i.ibb.co/j4XYzH9/Screenshot-1.png

Answer: a simple XSLT generating a shell script from the Sample names:

https://gist.github.com/lindenb/1d23b5d53f1a3694bbd1a7631935140f
Question: i am new the area of bioinformatics . Right now working on the tuxedo protocol . Trying to analyze a data based on the e coli reference genome.

this is my file location - **/Desktop/prasanth/TUXEDO/Experiments/Glucose_limited/alignersd/Bowtie/index$ ls
reference_genome.fa**

when i run the command -----**bowtie2 build ~/Desktop/prasanth/TUXEDO/Experiments/Glucose_limited/alignersd/Bowtie/index/reference_genome.fa**

This is the error that comes up

**Could not locate a Bowtie index corresponding to basename "build"
Error: Encountered internal Bowtie 2 exception (#1)
Command: /usr/bin/bowtie2-align --wrapper basic-0 build /home/ubuntu08/Desktop/prasanth/TUXEDO/Experiments/Glucose_limited/alignersd/Bowtie/index/reference_genome.fa
bowtie2-align exited with value 1**

i understand the problem with the location of file. But i am not able to solve the issue even after maintaining the reference genome in the same folder ..

ps: i am a beginner to the area of bioinformatics

Answer: Microfuge is right. Your command should be

bowtie2-build -f ~/Desktop/prasanth/TUXEDO/Experiments/Glucose_limited/alignersd/Bowtie/index/reference_genome.fa ~/Desktop/prasanth/TUXEDO/Experiments/Glucose_limited/alignersd/Bowtie/index/reference_genome

Then ~/Desktop/prasanth/TUXEDO/Experiments/Glucose_limited/alignersd/Bowtie/index/reference_genome is the address of your index. You can provide it at the time of alignment.
Question: Dear all,

please would you advise, based on your experience, which algorithm is more reliable in confidently calling INDELS in cancer genomes (I am referring especially to INDELS of length > 50bp, let's say 50bp-500bp long).

many thanks,

-- bogdan

Answer: pindel can do it but, when you get near to the 500bp mark, are you not then in the realm of small somatic copy number alterations? BBMap can also do this: https://www.biostars.org/p/254305/

In all cases, read length and insert size are obviously key.

Kevin
Question: Hi

I have a folder with 4000 files in it. I would like make a new folder, and copy into it all of the files which contain >50 fasta sequences. How do I do this?

I know that I need to create a simple loop, then grep '>' | wc -l and select only those with >50. But I am new to programming and am unsure how to write this properly.

Answer: Here is how

grep -c '^>' *.fasta | sed 's/:/ /' | awk '$2>50 {print $1}' | xargs mv -t new_folder/

**Update**

grep -c '^>' *.fasta | sed 's/:/ /' | awk '$2>50 {print $1}' | xargs cp new_folder/

P.S: Updated as @[genomax2][1] pointed out.

[1]: https://www.biostars.org/u/18713/
Question: Hi all,
here are some questions about paired-end sequencing for NGS:

What are the main differences between mate-paired sequencing and paired-end sequencing; Should I care when I use tools like 'samtools', maq, etc.... ? Should one, and only one short read, should be paired with another one (1-1)?

What is removing duplicates ? does it mean that a pair of short reads has been mapped at two distint positions on the genome or does it mean that a pair matched too many time at one position ?

Knowing that bwa sampe "Generates alignments in the SAM format given paired-end reads. Repetitive read pairs will be placed randomly", is there any need to "remove the duplicates" ?

How does picard MarkDuplicates work ? How can I find the reads that have been 'tagged' ? will it remove the reads from the BAM file ?

Thanks

Pierre

Answer: paired-ends and mate pairs are different protocols. The distance between mate pairs are much longer (2-5kb), while paired-end fragments are rarely more than 500bp apart and can even have negative distance (overlapping pairs)
Question: I want to query the molecular weight of multiple proteins from R given uniprot ids. I thought of using biomaRt or UniProt.ws, but didn't see which field I should query.

There are online tools performing the calculation given a uniprot ID (PIR, ExPASy), but you have to manually enter the input and save the output.

I want the MW that appears in UniProt when looking in "Sequences" under Mass (for example for the entry http://www.uniprot.org/uniprot/Q9Z2Y8)

Answer: You can get this information from the [TogoWS REST service][1], starting from UniProt ID. For example, using PROSC_MOUSE (which is the ID, Q9Z2Y8 is the *accession*):

```
library(rjson)
library(RCurl)

u <- getURL("http://togows.dbcls.jp/entry/uniprot/PROSC_MOUSE.json")
j <- fromJSON(u)

j[[1]]$sq$MW
# [1] 30049
```

Conversion from accession to ID should be easy using R/biomaRt (search Biostar for numerous examples).

[1]: http://togows.dbcls.jp/site/en/rest.html
Question: Hi,

I want to extract sequence with header from a fasta file with specific ID given in another file.

As for example:

File: `main.fa`

```
>ADC37925 pep:novel supercontig:GCA_000025145.2:CP001844:1841991:1848551:-1
MNLLKKNKYSIRKYKVGIFSTLIGTVLLLSNPNGAQALTTDNNVQSDTNQATPVNSQDTNVANNRGLANSAQNTPNQSATTNQGIFSTLIGTVLLLSNPNGAQALTTDNNVQSDTNQATPVNSQDTNVANNRGLANSAQNTPNQSATSNPNGAQAL
>EFB95474 pep:novel supercontig:GCA_000175495.1:cont1.9:99969:106529:-1
MNLLKKNKYSIRKYKVGIFSTLIGTVLLLSNPNGAQALTTDNNVQSDTNQATPVNSQDTNVANNRGLANSAQNTPNQSATTNQSTNQALVNHNNGSIANQGIFSTLIGTVLLLSNPNGAQALTTDNNVQSDTNQATPVGIFSTLIGTVLLLSNPNGAQ
>EFC04694 pep:novel supercontig:GCA_000175955.1:cont1.2:270427:276987:-1
MNLLKKNKYSIRKYKVGIFSTLIGTVLLLSNPNGAQALTTDNNVQSDTNQATPVNSQDTNVANNRGLANSAQNTPNQSATTNQSTNQALVNHNNGSIANQATPTSVQSSTPSAQNNNHTDGNTTATETVSNANNKDVVSNNTTLNVPNKTNENGS
>EFH37336 pep:novel supercontig:GCA_000178015.1:cont1.4:98713:105273:1
MNLLKKNKYSIRKYKVGIFSTLIGTVLLLSNPNGAQALTTDNNVQSDTNQATPVNSQDTNVANNRGLANSAQNTPNQSATTNQSTNQALVNHNNGSIANQVLLLSNPNGAQALTTDNNVQSDTNQATPVNSQDTVLLLSNPNGAQALTTDNNVQSD
```

File: `id.txt`

```
ADC37925
EFC04694
```

Expected Outcome:

```
>ADC37925 pep:novel supercontig:GCA_000025145.2:CP001844:1841991:1848551:-1
MNLLKKNKYSIRKYKVGIFSTLIGTVLLLSNPNGAQALTTDNNVQSDTNQATPVNSQDTNVANNRGLANSAQNTPNQSATTNQGIFSTLIGTVLLLSNPNGAQALTTDNNVQSDTNQATPVNSQDTNVANNRGLANSAQNTPNQSATSNPNGAQAL
>EFC04694 pep:novel supercontig:GCA_000175955.1:cont1.2:270427:276987:-1
MNLLKKNKYSIRKYKVGIFSTLIGTVLLLSNPNGAQALTTDNNVQSDTNQATPVNSQDTNVANNRGLANSAQNTPNQSATTNQSTNQALVNHNNGSIANQATPTSVQSSTPSAQNNNHTDGNTTATETVSNANNKDVVSNNTTLNVPNKTNENGS
```

Tried Perl Program :

I have found and made the Perl code work for a single input (ID) [[Source][1]]

```
########################## script.pl
#Usage: `perl script.pl file(s) 'searchTerm [searchTerm]' [>outFile]`

use strict;
use warnings;

my $term = join '.', map "\Q$_\E", split ' ', pop;
my $found;

while (<>) {
if (/^>/) {
$found = /$term/i ? 1 : 0;
}
print if $found;
}

###########################
```

Now how to use this perl code by taking the searchTerm from `id.txt` file to search in main.fa and give the outcome as above?

Any help will be very much appreciated.

I am new in Perl.

[1]: https://www.biostars.org/p/64149/

Answer: My Solution:

UCGS utilities

$ ./faSomeRecords main.fasta id.txt output.fa

option `-exclude` will output sequences not present in "main.fasta"
Question: Is there any way to detect if your rna-seq data is unstranded or stranded

Answer: In case anybody currently looking comes across this post...

The easy-to-use Salmon will check for you as described [here][1]. You can see what the resulting abbreviations correspond to with a nice illustration [here][2].

[1]: https://salmon.readthedocs.io/en/latest/salmon.html#what-s-this-libtype
[2]: https://salmon.readthedocs.io/en/latest/library_type.html#fraglibtype
Question: Hi folks. I need to run a de novo short-read genome assembler (on a paired-end/mate-pair library) that prefers outputting shorter but error-free contigs rather than longer contigs/scaffolds which may be mis-assembled. What assembler or what specific setting in an assembler of choice do you recommend to yield such contigs (as error-free as possible and no contig overlappings)?

Answer: According to the first GAGE paper, SGA makes shorter, but very much correct contigs. See http://genome.cshlp.org/content/early/2012/01/12/gr.131383.111.full.pdf