cannin/biostars_qa huggingface dataset:
Question: Hi all,

I used bamToFastq to convert my bamfiles into R1 and R2 fastq files. I tried samtools fastq before but found it was erroneous (so many missing sequences). However, I still have an issue with bamTofastq from bedtools. When I add the read counts obtained from ```grep -c "@" R1.fastq``` and ```grep -c "@" R2.fastq```, it is always slightly less than the count from the bamfile ```samtools view -c in.bamfile```. Why might this be the case? I haven't found any documentation to suggest that this is a normal thing. The R1 and R2 fastq counts should equal the counts in the bamfile so what am I doing wrong with the conversion??

Thank you.

Answer: > it is always slightly less than the count from the bamfile samtools view -c in.bamfile. Why might this be the case?

you bam contains supplementary, secondary alignments.
Question: Hi all,

I'm facing a very annoying error in R while assigning row names to my data matrix. I have some RNA-seq data that I'm considering clustering in R. I'm using gene names as row names for my expression matrix but it keeps reporting that there are duplicate names. Some un-annotated genes have been assigned with some IDs that start with numbers. I don't understand how to deal with this error? Is there a way to work around it? because I cant change the gene names.

EDIT:

gene sample1 sample2 sample3
Mar-01 4.19504 3.9006 4.15683
Mar-02 3.0554 3.4261 3.76675
un_A_2 1.1515 1.2455 0.563484
un_A_3 98.2504 120.341 101.753
ENSGALG00000008227 39.6383 12.8651 38.2281
ENSGALG00000008242 5.71557 7.79314 9.40917
ENSGALG00000008277 24.6231 28.3207 24.9288
CNN3 141.708 134.476 144.514
CNNM1 0.840218 0.963683 0.619086
CNNM2 16.0282 12.1301 12.4665

Many thanks.

Answer: One way of dealing with this is in R is the function make.names with the option unique=TRUE, see ?make.names.

> nams = c("bl-a","bl-a","bl-a", "foo" )
> df = data.frame(matrix (1:4))
> df
matrix.1.4.
1 1
2 2
3 3
4 4
> rownames(df) = nams
Error in `row.names<-.data.frame`(`*tmp*`, value = value) :
duplicate 'row.names' are not allowed
In addition: Warning message:
non-unique value when setting 'row.names': ‘bl-a’
> rownames(df) = make.names(nams, unique=TRUE)
> df
matrix.1.4.
bl.a 1
bl.a.1 2
bl.a.2 3
foo 4
Question: I'd like to extract consensus of reads from a bam/sam file which would follow particular criteria:

1. The portions of reference which have 0 coverage are represented in
the consensus as N
2. If there is at least single read mapping to a particular nucleotide,
this nucleotide is retrieved in a consensus.
3. If there is conflicting evidence from reads, the nucleotide with
highest frequency (lets say 75%) is retrieved in consensus. If there is no nucleotide with >= 75% frequency, N is retrieved in the consensus.

Is there a tool which would enable this or similar level of fine-tuning of the consensus extraction from bam/sam?

What I tried: I used samtools, bcftools and vcfutils to get consensus:

samtools mpileup -uf reference.nt mapping.bam | bcftools call -c | vcfutils.pl vcf2fq > consensus.fasta

However, the consensus I'm getting this way contains a substantial proportion of the reference sequence and I can't find parametres which could modify how the consensus is calculated.

Yet another approach - extracting consensus using IGV viewer (option Copy consensus sequence) gives me basically the consensus I'm looking for but there is apparently no way how to automatize this for hundreds of reference sequences which I have.

Note: I posted this question also as a comment to my original question [Map eukaryotic genomic reads to transcript reference of closely related organism][1] however I'm posting it here separately as it diverged from my original question.

[1]: https://www.biostars.org/p/288924/#289983

Answer: simple-consensus-per-read-group by Thomas Sibley is the tool I was looking for (https://github.com/MullinsLab/simple-consensus-per-read-group). It was recently updated so it can now correctly process BAMs with multiple scaffolds/contigs.
Question: How can I construct a phylogenetic tree based on the SNP's shared between strains? I have whole genome SNP calls for 10 different strains in a multi-sample vcf.

Are there any tools that can take the vcf as an input for creating phylogenetic trees? Or do I need to convert the multi-sample vcf to another matrix? Which kind of matrix would that be how can I create it from the vcf?

Is there a list somewhere of popular packages that kan be used for creating phylogenetic trees? Or guides on how to go from a multi-sample vcf to a plylogenetic tree.

Answer: Here is what I did in the SNPrelate package to get a dendogram and pca from my multisample vcf file

#vcf to GDS
snpgdsVCF2GDS("my.vcf", "my.gds")
snpgdsSummary("my.gds")
genofile <- openfn.gds("my.gds")

#dendogram
dissMatrix <- snpgdsDiss(genofile , sample.id=NULL, snp.id=NULL, autosome.only=TRUE,remove.monosnp=TRUE, maf=NaN, missing.rate=NaN, num.thread=10, verbose=TRUE)
snpHCluster <- snpgdsHCluster(dist, sample.id=NULL, need.mat=TRUE, hang=0.25)
cutTree <- snpgdsCutTree(snpHCluster, z.threshold=15, outlier.n=5, n.perm = 5000, samp.group=NULL,col.outlier="red", col.list=NULL, pch.outlier=4, pch.list=NULL,label.H=FALSE, label.Z=TRUE, verbose=TRUE)

#pca
sample.id <- read.gdsn(index.gdsn(genofile, "sample.id"))
pop_code <- read.gdsn(index.gdsn(genofile, "sample.id")
pca <- snpgdsPCA(genofile)
tab <- data.frame(sample.id = pca$sample.id,pop = factor(pop_code)[match(pca$sample.id, sample.id)],EV1 = pca$eigenvect[,1],EV2 = pca$eigenvect[,2],stringsAsFactors = FALSE)
plot(tab$EV2, tab$EV1, col=as.integer(tab$pop),xlab="eigenvector 2", ylab="eigenvector 1")
legend("topleft", legend=levels(tab$pop), pch="o", col=1:nlevels(tab$pop))
Question: Does anyone know of any software package that can convert .impute2 data to .mldose (i.e. imputed data from IMPUTE2 to imputed data from MACH)? I have tried impute2mach in GenABEL, but it has repeatedly failed with a known error...

Direct conversion preferred, but I'm happy with indirect as long as it works!

Answer: I have written a brief cookbook to perform this conversion in UNIX:

http://openwetware.org/wiki/User:Jonathan_R._I._Coleman/Notebook/Notes_and_Protocols/2014/06/27
Question: Hello everybody,

Here i my issue, I have a Trinity.fasta file like this

>TRINITY_DN5631_c0_g2_i1 len=947 path=[0:0-946]
TACAACTTGAACATCAACAATGGTTGCGCAGCTATTGCCATCCGCGACGTTCGAGGACTGCGTGCGAA
>TRINITY_DN62279_c1_g1_i1 len=298 path=[0:0-297]
TATTACCATTATTATTATTATCATATTTATGTTCATTATTATCATTATCATAATCATTATCATCTTGATA
...
And I also have a list of id in a id.txt file :

TRINITY_DN16359_c0_g1_i4
TRINITY_DN62279_c1_g1_i1
...

I am trying to extract for my fasta file the sequences that have their id in my txt.file.
I am using seqkit for that, but with no success :

seqkit grep -n -f id.txt Trinity.fasta -o result.fa

Does anyone know how to fix it ?

Answer: Referring to the [usage](https://bioinf.shenwei.me/seqkit/usage/#grep), you should not switch on `-n` , just use `seqkit grep -f id.txt seqs.fa`.

-n, --by-name match by full name instead of just id
-i, --ignore-case ignore case
-r, --use-regexp patterns are regular expression

For `-nrif`, it's **partly matching** full FASTA header by regular expression and case-ignored, with patterns from file.
*This may produce some unwanted results*. For example, `seq_1` matches `seq_10` with `-nri`.
Question: In the documentation of Gemini there is that you should

1. [Decompose][1] the original VCF such that variants with multiple alleles are expanded into distinct variant records; one record for each REF/ALT combination.
2. [Normalize][2] the decomposed VCF so that variants are left aligned and represented using the most parsimonious alleles.

http://gemini.readthedocs.org/en/latest/index.html

This sound like a good thing to do because is makes it easier to asign correct IDs and to compare variants.

Only problem I have seen mentioned is that for samples that have an ALT1, ALT2 genotype your genotype is now split over 2 vcf records; MISSING, ALT1 and MISSING , ALT2 . Or even REF, ALT1 and REF, ALT2 . Both don't correctly represent the genotype of the sample.

[1]: http://genome.sph.umich.edu/wiki/Vt#Decompose
[2]: http://genome.sph.umich.edu/wiki/Vt#Normalization

Answer: There are a few ways to skin this cat, and it is also an area with fairly active development. The central difficulty is that there often are multiple ways to represent the same variant in VCF, particularly in cases where block substitutions or indels are involved, and there is no "right" representation.

The decomposition/normalization approach has the downside that the process tends to destroy a lot of the good information that is contained in the original call set (e.g. phasing information, INFO/FORMAT annotations, quality scores). In addition, even after decomposition the results can be arbitrary (and so may not match up with with the coordinates you are getting your IDs from anyway, defeating the purpose).

An alternative approach is to have smarter comparison tools which are directly aware of representational ambiguity, by performing variant comparison at the haplotype level. AFAIK [CGI calldiff][1] and [RTG vcfeval][2] were independently the first to implement this strategy, and new tools are finally catching on, in varying stages of development ([SMaSH][3], [vgraph][4], [hap.py][5]). These tools replay the variants from the VCF into the reference and determine whether variants match by whether the resulting haplotypes match. With vcfeval the full VCF annotation information is preserved during the comparison (not so with hap.py, vgraph doesn't currently output VCF, and I haven't used calldiff or SMaSH)

In particular, the haplotype comparison tools are the current state of the art for same-sample call-set comparison (either between callers, or comparing with a benchmark set) -- certainly in the case of vcfeval this was the motivating driver in the development. The decomposition/normalization approach is more useful if you want to establish a population-level database where variants are converted to a "canonical" form with limited annotation requirements. Of course there is nothing to say you cannot use both techniques, depending on what you are trying to achieve.

[1]: http://www.completegenomics.com/public-data/analysis-tools/cgatools/
[2]: http://realtimegenomics.com/products/rtg-tools/
[3]: http://smash.cs.berkeley.edu/
[4]: https://github.com/bioinformed/vgraph
[5]: https://github.com/Illumina/hap.py
Question: Hello,

does anyone know the code for each color of UCSC genome browser? Is there a list or something?

For example this one "color=250,0,0" gives me the color RED and this one "color=0,0,250" is blue.

Is there a list in which i can select the serial number for other colors?

Thanks in advance

Answer: It uses standard [RGB codes][1].

You can Google them to get the code for the color you want (here is first hit: https://www.rapidtables.com/web/color/RGB_Color.html)

[1]: https://en.wikipedia.org/wiki/RGB_color_model
Question: The mouse gene: ENSMUSG00000004500 has 2 human homologs: ENSG00000249471, and ENSG00000083812. Looking the opposite direction, the indicated mouse gene is the only mouse homolog of these two human genes. I was hoping someone from Ensembl could clarify why is this labeled "many2many" instead of "one2many"?

I retrieved a table of mouse homologs for human genes (Ensembl v104, data set Human Genes GRCh39) using the following biomart query:

Attribute set: "homologs radio button"
Attributes:

1. Gene.stable.ID
2. Mouse.gene.stable.ID
3. Mouse.homology.type
4. Homology related stats [ %identities, GOC, Alignment Coverage, and Confidence]

As I was reviewing the table, I found several records where the Mouse.homology.type was labeled "ortholog_many2many" where the homologous mouse gene had multiple human homologs, but those human homologs are only associated with the single mouse gene. Shouldn't these genes be labeled "one2many" as is true for most other similar cases:

R Filtering Code:

> with(
+ hom_2_mus %>%
+ select(
+ Human_ID=Gene.stable.ID,
+ Mouse_ID=Mouse.gene.stable.ID,
+ Mouse.homology.type
+ ) %>%
+ group_by(Human_ID) %>%
+ filter(n() == 1) %>% # Human genes with only one mouse homolog
+ group_by(Mouse_ID) %>%
+ filter(n() > 1), # Mouse genes with multiple human homologs
+ table(Mouse.homology.type)
+ )
Mouse.homology.type
ortholog_many2many ortholog_one2many
44 1477

Answer: Hi Adam,

This particular mouse gene is the only homologue of the two human gene stable IDs, but has a listing of "many2many" still because this mouse gene has undergone a duplication event following the speciation. Although the homologous relationship is only shown between the gene stable IDs - this mouse gene has more than one protein coding gene attached, and so effectively the two human genes show a "many2many" relationship on the gene tree - with the 2 mouse proteins. It could be noted that the human side homologues include more than just those two genes, and those two gene stable IDs actually encode multiple protein coding genes.
Question: Hi - I'm new to Plink and am trying to read in a transposed fileset so that I can then convert it to bed/bim/fam, but I keep getting this message that says "Note: Variant ####### is triallelic. Setting rarest alleles to missing." Where there are many lines and the ####### ranges from 0-497197. But it still creates files that have the extension `.temporary.bed.tmp`, `.temporary.bim`, and `.temporary.fam`.

So my problem is that I can't figure out why it thinks that they're all triallelic when I can look at my .tfam file and see that they aren't (or at least don't appear so to me). Does anyone have any suggestions?

My tped input looks like this (just with 497198 SNPs):

```
12 rs1000000 0 126890980 A A B B B B B B
4 rs10000023 0 95733906 B B B A B A B A
4 rs10000030 0 103374154 B B B B B B A B
4 rs10000041 0 165621955 A A B B B B B A
4 rs10000042 0 5237152 B B B B B B B B
```

My plink command looks like this:

./plink --tfile myfile --recode --out myfile

And the return I get from plink looks like this:

Note: Variant XXXXXX is triallelic. Setting rarest alleles to missing.

In addition, I'm also getting the error (on the very last line) but I thought the centimorgan position was allowed to be 0?:

Error: Invalid centimorgan position on line 2 of .tped file

Any help would be much appreciated - thanks in advance!

Answer: The `.tped` loader did not properly detect when a line didn't have enough genotypes for the number of samples in the `.tfam` file. This error is properly reported in the Nov 23 development build.

You now need to check why your `.tfam` has more lines than there are genotype pairs in the `.tped`.
Question: I have a large set of footprint intervals that range from 11 to 25bp For the purpose of motif discovery I would like to extend all intervals to, for example, 50bp. Intervals should be extended equally from both sides. I would usually use 'bedtools slop' for fixed length intervals, but this would not appear to work with variable length.

It would be great if anyone could advise me how to use bedtools, or something else. I have a nagging feeling I am missing something obvious, so apologies in advance!

Answer: Here's a way that I think should extend both ends of BED elements to the desired target length:

$ TARGET_LENGTH=50
$ awk -vF=${TARGET_LENGTH} 'BEGIN{ OFS="\t"; }{ len=$3-$2; diff=F-len; flank=int(diff/2); upflank=downflank=flank; if (diff%2==1) { downflank++; }; print $1, $2-upflank, $3+downflank; }' in.bed | sort-bed - > out.bed

Non-even length elements or a non-even target length will require flank lengths that are unequal. Sounds like this is not a problem.

You might adjust the logic to randomly pick which of `upflank` or `downflank` to decrement or increment in this case, so that you don't impart a bias from this adjustment (esp. if original elements are stranded, like footprints that will ultimately be mapped to TF binding sites or other stranded elements), e.g.:

$ TARGET_LENGTH=50
$ awk -vF=${TARGET_LENGTH} 'BEGIN{ OFS="\t"; }{ len=$3-$2; diff=F-len; flank=int(diff/2); upflank=downflank=flank; if (diff%2==1) { if (rand() >= 0.5) { downflank++; } else { upflank--; } }; print $1, $2-upflank, $3+downflank; }' in.bed | sort-bed - > out.bed
Question: Hello,
I have a fasta file with lines of format -

>FBti0019256 type=transposable_element; loc=2L:22300300..22304444; name=invader2{}555; dbxref=FlyBase_Annotation_IDs:TE19256,FlyBase:FBti0019256; MD5=d9259a0e33aad699215e64916bd47a5b; length=4145; release=r6.19; species=Dmel;

I would like to convert these lines into a bed file of format -

chr2L /t 22300300 /t 22304444 /t invader2

Is there a program that can directly perform this conversion or an awk command that can do this easier? Please let me know, thank you for your help.

Answer: Assumption is that chromosomes are numbered, no x and y chromosomes. For fasta headers only, user can use grep

$ sed 's/.*loc.\(\b.*\b\):\(\b.*\b\)\.\.\(\b.*\b\);.*=\(\b.*\b\){.*/chr\1\t\2\t\3\t\4/g' test.txt
chr2L	22300300	22304444	invader2

$ cut -f2,3 -d";" test.txt| cut -d= -f2,3 | awk -v OFS="\t" -F':|=|;|\\..|{' '{print "chr"$1,$2,$3,$5}'
chr2L	22300300	22304444	invader2

$ cat test.txt
>FBti0019256 type=transposable_element; loc=2L:22300300..22304444; name=invader2{}555; dbxref=FlyBase_Annotation_IDs:TE19256,FlyBase:FBti0019256; MD5=d9259a0e33aad699215e64916bd47a5b; length=4145; release=r6.19; species=Dmel;

1000 bash cuts:

$ cut -f2,3 -d";" test.txt| cut -d= -f2,3 | cut -f1,2 -d: --output-delimiter=$'\t'| cut -f1,2,3 -d'.' --output-delimiter=$'\t' | cut -f1,2 -d';' --output-delimiter=$'\t'| cut -f1,2 -d'=' --output-delimiter=$'\t' | cut -f1,2 -d'{' --output-delimiter=$'\t'| cut -f1,2,3,4,6

2L	22300300	22304444	invader2
Question: Hi,

I'm trying to extract reads based on their start coordinate in a bam file. I've tried using samtools view but that seems to give all reads covering that region, not originating there.

Apologies if this has been asked elsewhere but I couldn't find an answer via Google.

Many thanks.

Answer: Something like:

`samtools view -h -q 10 in.bam chr1:1000-1500 | awk 'BEGIN{OFS="\t"}{if($1 ~ @) {print} else {if($4 >= 1000 || $4 <=1500) {print} else {}}}' | samtools view -Sbo output.bam -`

One could restructure this to get around the default print action in awk, but this will also get the job done.
Question: I have a set of vcf files that were filtered using GATK hard filtering. I filtered the snps and then the indels seperately and then merged the two and got a vcf with a list of polymorphisms in which the polymorphisms that had failed the filters weere marked as such (of course) Now I would like to make a vcf that lacks the snps and Indels that failed the filtering. What command should I run.

Answer: If you're referring to the `PASS` flag, you can use anything from `vcftools` to plain `awk` to plainer `grep`. If there's more to the PASS criteria than just the flag, you're going to need to elaborate on that.

Also, do you mean `variants` when you say `polymorphisms`? I know people may use the terms interchangeably, but they do not mean the same thing. `Variants` are the most generic descriptor; they refer to all loci where multiple alleles are found, whereas `polymorphisms` assume the functional impact of the variant to result in a polymorphic phenotype that is not usually pathogenic.
Question: hi everyone,

I have just started learning genomics as a part of my bioinformatics degree and I've been introduced to linux for handling fastq files and using fast qc. can anybody suggest some good learning resources which is more inclined towards linux for genomics.

Answer: My advice would be (to echo Ram) focus on a subset of small command line commands and don't worry too much about learning about the nuts and bolts of Linux. In a rough order of importance I would say learn the following:

1. ``cd``, ``cat``, ``mkdir``, ``ls -lh``, ``rm``, ``touch``, ``head``, ``tail``, ``echo``, ``find``, ``man``, ``../``, ``sudo``
2. ``cut``, ``for``, ``while``, ``if else``
3. ``grep``, ``sed``, ``awk``, ``parallel``
Question: Hi all,

I have a couple of fastq files containing reads starting with different name like:
@HWI-ST865:463:C7C8KACXX:2:2316:21016:100943 1:N:0:TAAGGCGA
@HWI-ST1178:227:C7C95ACXX:7:1101:1581:2125 1:N:0:TAAGGCGA

My question is: how can I split them in two parts?
I tried to use some tools like fastx_toolkit but I cannot create a proper barcode file
Is there any easy way to do that such as a grep command, cause i also tried with grep but i got an output containing only the first line of the reads and missed the other three

Thank you in advance!

Answer: You can use either Heng Li's `bioawk` or `grep -A 3`. The former is a wrapper on awk to make it work with separators used in biological data formats, and the latter is a grep that picks up the matching line+3 lines that follow.
Question: I'm looking to use something different than the built-in color palettes in `Seurat`, but I'm finding it very difficult to find packages that can handle 20+ discrete colors. I've looked through `RColorBrewer` and `colorspace`, as well as `ggsci`, and all max out around 12 colors from what I can see. Does anyone know of other color packages that can handle a higher number of colors, or am I best off using `scale_color_manual()` and doing it myself in `ggplot2`?

Answer: > am I best off using scale_color_manual() and doing it myself in ggplot2

Yes.

Keep in mind though that there is a reason for the limit of discrete colors! It's not the computers or tools that have a problem distinguishing more than 12 discrete colors, it's the human perception that tends to be the limiting factor here!

That being said -- I don't know why you find that `RColorBrewer` is limited.

```
> library(RColorBrewer)
> colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(20)
[1] "#4575B4" "#5D8CC0" "#75A3CC" "#8DBBD8" "#A5CCE2" "#BEDDEB" "#D7EDF4" "#E6F5EC" "#F0F9DA" "#FAFDC8" "#FEFAB7"
[12] "#FEF0A8" "#FEE699" "#FDD78A" "#FDBD78" "#FCA267" "#FA8856" "#EE6A46" "#E24D36" "#D73027"
```
Question: I would like a list of drug-variant interactions (i.e. a variant for which its mutational status affects the efficacy of a particular drug).

I noticed that [nightly-ClinicalEvidenceSummaries.tsv][1] currently contains 593 drug-variant interactions (as of March 22, 2016). Is this the full list of all drug-variant interactions currently contained in the CIViC database?

If so, great! Otherwise, it seems like the alternative is to use the API to obtain variant/evidence information and then parse the resulting JSON. Thanks.

[1]: https://civic.genome.wustl.edu/#/releases

Answer: That sound right. The nightly TSV file should be a complete representation of the total number of evidence statements in CIViC. However, the level of detail is reduced compared to what you could get from the API. Also, I see that there are some formatting issues still occurring with the TSV file. We are working to fix this. However, the API is definitely recommended.

Check out this set of python tools (under active development) that pulls CIViC data via the API in various ways.
https://github.com/griffithlab/civic-api-client