cannin/biostars_qa huggingface dataset:
Question: This has been bugging me for awhile now. I know there are multiple mapping reads in my alignments, but bowtie2 does not assign the 256 flag for any of my reads. Bowtie2 says in the manual that is should be doing this:

> Each reported read or pair alignment beyond the first has the SAM
> ‘secondary’ bit (which equals 256) set in its FLAGS field.

But when I check for any secondary alignments:

$ samtools view -f 256 my_file.bam | wc -l
0

There are apparently no secondary alignments.

**How do I know there actually are multiple mapping reads?** Well bowtie2 also uses the `XS:i:` field

> XS:i: Alignment score for the best-scoring alignment found other
> than the alignment reported. Can be negative. Can be greater than 0 in
> --local mode (but not in --end-to-end mode). Only present if the SAM record is for an aligned read and more than one alignment was found
> for the read. Note that, when the read is part of a
> concordantly-aligned pair, this score could be greater than AS:i.

So, since `XS:i` is ***only*** present if more than one alignment was found for the read, let's check that:

$ samtools view my_file.bam | grep "XS:i:" | wc -l
42821156

Am I missing something here? There are clearly multiple alignments for a LOT of my reads, so why is there not a single alignment with 256 flag set?

Answer: Oops... I know what is happening now. Bowtie2, by default, only reports the best alignment. This is why I have zero secondary alignments in my bam file, despite having many "multi-mappers".

This information was in the **Reporting** section of the [bowtie2 manual][1]. I think the manual could be more clear... But I admit it's my fault for missing this before. In case anyone cares, you can use the `-k` option to report a given number of secondary alignments, or the `-a` option to report all alignments.

[1]: http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml#reporting
Question: Hello all,

The Broad Institute did the CNV analysis for the TCGA samples. The only information I have been able to find about the reference genome* does not give the version of hg19 they are using. I have two questions:

1. Is hg19 exactly the same as Ensembl's GRCh37?
2. It looks from the Ensembl site that the data has had various patches and updates. I have a list of CNV start end co-ordinates and want to extract the genes involved. How much difference will it make if I use a version with the wrong patches?

Thanks,
Stephanie

https://confluence.broadinstitute.org/display/GDAC/FAQ#FAQ-Q Whatreferencegenomebuildareyouusing

Answer: For the autosomes and sex chromosomes, hg19 and GRCh37 are the same. For the purposes of start and end coordinates, the patches won't make a difference for the standard chromosomes.
Question: I worked in DNA K-mers counting and I prepare this formulation to solve the counting using perfect hash table :
![enter image description here][1]
![enter image description here][2]

Using Rcpp API(C++) to integrate the code into R:

#include
using namespace Rcpp;
/*
this code can be used with c++ by replacing IntegerVector by std::vector
*/
//************************************************
inline const short int V (const char x){
switch(x){
case 'A':case 'a':
return 0;
break;
case 'C':case 'c':
return 1;
break;
case 'G':case 'g':
return 2;
break;
default:
return 3;
break;
}

}

inline unsigned int X0( const std::string A,const int k ,const int n){
unsigned int result=0;
int j=k;
for( int i=n-1;i>n-k-1;i--) {
result+= pow(4,k-j)*V(A[i]);
j--;
}
return result;
}

// [[Rcpp::export]]
inline IntegerVector kmer4(const std::string A,const int n,const int k)
{

IntegerVector P(pow(4,k));
int x=X0(A,k,n);
P[x]++;
const int N=pow(4,k-1);
for( int i=n-k-1;i>-1;i--){
x=N*V(A[i])+x/4-x%4/4;
P[x]++;
}
return P;
}

There are two questions:

Assuming the index x of a kmer,the compliment of x is then (4^k)-x-1. Can we get the reverse using a numeric operation like preceding formula?
There are two problems in running time: iteration over the string and vector creation where the k is over then 8. Are there ideas to solve these problems?

[1]: https://s15.postimg.org/dwqyrpbvv/Capture.jpg
[2]: https://s15.postimg.org/o7iozlzfv/Capture2.jpg

Answer: I think that the problem is solved using this code using incrementation.

inline unsigned int X1( const std::string A,const int k ){
unsigned int result=V(A[0]);
for( int i=1;i
Question: Hi all,

I have RNA-Seq time-series data generated by HTseq counts. There are sort of 2 replicates. The First replicate has 4 time point samples(day0,2,5,14), while the Second has 5 time point samples(day0,2,5,15,30). The data was loaded as a read.table matrix in R.
The following is the script used:

setwd("Downloads")
df_A <- read.table(file='Replicate_sample_Raw_RPM_counts-2',header=T,sep='\t')
df_B <- read.table(file='HT-seq_counts.txt-2',header=T,sep='\t')
merged_df <- merge(df_A,df_B,by='geneID')
write.table(merged_df,file='merged_count.txt',row.names=F,quote=F,sep='\t')
counts<-read.csv("merged_count.txt",header=T,sep="\t")
data<-counts[-1]
rlog(data)

head(data)

day0.x day2.x day5.x day14 day0.y day2.y day5.y day15 day30
1 358 422 241 617 145 508 389 357 594
2 11 31 44 26 8 24 41 49 49
3 7 3 33 392 2 5 25 159 155
4 26 45 74 5624 45 175 94 4604 14238
5 4 10 66 338 19 13 70 229 242
6 477 138 64 21 747 507 98 25 22

But when I try to to run rlog or rlogTransformation I get the following error:

Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘sizeFactors’ for signature ‘"data.frame"’

How can I improve the script to make rlog run so that I can have my HT-seq counts normalized?

Answer: I would highly recommend you read the [DESeq2 vignette][1].

You first need to create a design matrix. This will be a 3 columned table, where each row represents a replicate, and each column represents, time, donor/sample, and treatment (If you have control replicates. If all samples are treated with the same stimulation (i.e. no controls), then you don't need a column for it.).

You can proceed to create a DESeq data object (your data is currently a "raw matrix") from your count table and your design table.

dds<-DESeqDataSetFromMatrix(data, design, ~ time + sample + treatment)

Your design formula needs to take into account your time points as well as the two replicated samples.

Then, estimate the size factors (partial cause of the error message), and proceed to rlog the data-set

norm<-estimateSizeFactors(dds)
expr<-rlog(norm)

Also, you might want to look into VST normalizing your data (run it with the blind parameter = FALSE). VST normalization takes into account your experimental design when normalizing. You might want to remove the 30 day time point replicate, this will make your time-course data-set more even.

[1]: https://www.bioconductor.org/packages/3.3/bioc/vignettes/DESeq2/inst/doc/DESeq2.pdf
Question: I have paired end data in a SAM file and I want to extract all the reads which correspond to specific contig ids. My contig ids are in a text file with one ID per line.

Output of SAM/BAM is preferred.

Thanks for your help.

Answer: ah, I see. Using awk, matching chromosomes RF01 or RF02

/samtools view -h S1.bam | awk -F '\t' 'function fun(C) { return C=="RF01" || C=="RF02";} /^@/ {print;next;} {if(fun($3) || fun($7)) print;}'

or using samjdk: http://lindenb.github.io/jvarkit/SamJdk.html

java -jar dist/samjdk.jar -e 'Predicate f=C->C.equals("RF01") || C.equals("RF02"); return (!record.getReadUnmappedFlag() && f.test(record.getReferenceName())) || (!record.getMateUnmappedFlag() && f.test(record.getMateReferenceName()));' in.bam
Question: Dear all,
My question is regarding the making of an image that will give information about the position of domains on proteins.
I have thousands of proteins and for these proteins, we have got the domain type and their position after manual curation from many databases.
With this information, I want to make an image that will have all the domains of specific protein with their positions plotted on an image using any programming language viz-R, Python, HTML, etc. Any help in this regards will be highly appreciated.

Protein Start End Domain
**Aco021107	1	1000	Full protein length**
Aco021107	496	702	DUF1336
Aco021107	169	337	START
Aco021107	14	112	PH

Thanks in advance.

Answer: Haven't tried this, but it seems to be exactly what you want:

https://f1000research.com/articles/7-1105/v1

A very specific solution you may need to fiddle with to suit your needs:

http://rforbiochemists.blogspot.com/2015/11/drawing-protein-domain-structure-using-r.html

The links below do not offer a programmatic solution, but they still may be helpful.

Pfam has its own utility based on JSON to create graphics:

https://pfam.xfam.org/generate_graphic

Other choices:

- https://prosite.expasy.org/mydomains
- http://ibs.biocuckoo.org/
Question: Hi everyone,

I am bioinformatician and I worked on "Mathematical Modeling and Simulation" of signaling pathways. Now, I want to work on Microarray data. I am new in this field, and I didn't know more about microarray data. How we generate (criteria) the excel file, which includes (Gene name,Probe.., Intensity, log ratio, etc) from our dataset. If the data is already normalized then what will be the next step? **Is there any good tutorial for beginner, which guides step by step from data selection to analysis**. I need the valuable suggestions from all the experts in this field.

Thank you

Answer: If you are beginner in microarrays, please follow the limma user's guide.

http://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf
Question: After spending numerous hours trying to install bcl2fastq 2.20 in my server I decided to ask for help.

I am trying to install bcl2fastq 2.20 in Ubuntu 18.04. I have followed multiple recommendations for fixing inconveniences while installing bcl2fastq with newer boost versions (as in these posts [here][1] and [here][2] ).

I don't get any errors after running configure from the build directory. My error appears after typing make afterwards.

$ make
Scanning dependencies of target BCL2FASTQ_OPT
[ 0%] Built target BCL2FASTQ_OPT
Scanning dependencies of target bcl2fastq_common
[ 1%] Building CXX object cxx/lib/common/CMakeFiles/bcl2fastq_common.dir/Logger.cpp.o
/media/software/bcl2fastq/bcl2fastq/src/cxx/lib/common/Logger.cpp: In constructor ‘bcl2fastq::common::detail::LogStream::LogStream()’:
/media/software/bcl2fastq/bcl2fastq/src/cxx/lib/common/Logger.cpp:122:13: error: ‘clog’ is not a member of ‘std’
, ias_(std::clog)
^~~~
/media/software/bcl2fastq/bcl2fastq/src/cxx/lib/common/Logger.cpp:122:13: note: suggested alternative: ‘log’
, ias_(std::clog)
^~~~
log
/media/software/bcl2fastq/bcl2fastq/src/cxx/lib/common/Logger.cpp: In constructor ‘bcl2fastq::common::detail::LogStream::LogStream(bcl2fastq::common::LogLevel::value_type)’:
/media/software/bcl2fastq/bcl2fastq/src/cxx/lib/common/Logger.cpp:130:13: error: ‘clog’ is not a member of ‘std’
, ias_(std::clog)
^~~~
/media/software/bcl2fastq/bcl2fastq/src/cxx/lib/common/Logger.cpp:130:13: note: suggested alternative: ‘log’
, ias_(std::clog)
^~~~
log
/media/software/bcl2fastq/bcl2fastq/src/cxx/lib/common/Logger.cpp: In copy constructor ‘bcl2fastq::common::detail::LogStream::LogStream(const bcl2fastq::common::detail::LogStream&)’:
/media/software/bcl2fastq/bcl2fastq/src/cxx/lib/common/Logger.cpp:146:13: error: ‘clog’ is not a member of ‘std’
, ias_(std::clog)
^~~~
/media/software/bcl2fastq/bcl2fastq/src/cxx/lib/common/Logger.cpp:146:13: note: suggested alternative: ‘log’
, ias_(std::clog)
^~~~
log
cxx/lib/common/CMakeFiles/bcl2fastq_common.dir/build.make:57: recipe for target 'cxx/lib/common/CMakeFiles/bcl2fastq_common.dir/Logger.cpp.o' failed
make[2]: *** [cxx/lib/common/CMakeFiles/bcl2fastq_common.dir/Logger.cpp.o] Error 1
CMakeFiles/Makefile2:257: recipe for target 'cxx/lib/common/CMakeFiles/bcl2fastq_common.dir/all' failed
make[1]: *** [cxx/lib/common/CMakeFiles/bcl2fastq_common.dir/all] Error 2
Makefile:126: recipe for target 'all' failed
make: *** [all] Error 2

ANy ideas on how to fix this?
Suggestions are highly appreciated.

[1]: https://backwardincompatible.com/post/169360794395/compiling-illumina-bcl2fastq-220-on-ubuntu-with
[2]: https://gist.github.com/jblachly/f8dc0f328d66659d9ee005548a5a2d2e

Answer: I finally found a post on how to install bcl2fastq v2.20 from a tarball. I followed the instructions and installed it with no errors!!! https://sarahpenir.github.io/linux/Installing-bcl2fastq/
Question: I have a multi-sample vcf file and I want to get a table of IDs on the left column with the variants in which they have an alternate allele in. It should look like this:

> ID1 chr2:87432:A:T_0/1 chr10:43234:C:G_1/1
ID2 chr2:87432_A:T_1/1
> ID3 chr11:432434:T:G chr14:34234234:C:G chr20:34324234:T:C

This is to then read into R

I have tried combinations of:

bcftools query -f '[%SAMPLE\t] %CHROM:%POS:%REF:%ALT[%GT]\n'

but I keep getting sample IDs overlapping on the same line and I can't quite figure out the sytnax.

Your help would be much appreciated

Answer: something like a loop over each sample:

```
bcftools query -l in.vcf.gz | while read S; do echo -n "$S" && bcftools view --samples $S -O u in.vcf.gz | bcftools query -i 'AC>0' -f '\t%CHROM:%POS:[%GT]' && echo ; done
```

I'm too lazy to convert the alleles' indexes to their value...
Question: I have been trying to find if TCGA PRAD (prostate cancer) RNA-seq data strand specific or not. I can not find any relevant info. I did find that it is Illumina protocol, but nothing more. Does anyone if they are strand specific? Do you have the right link?
Thanks

Answer: The short answer is no.

The specific RNA-Seq library preparation kit TCGA used actually "could" support stranded library; however, it was used in a manner that libraries are actually strandless.
Question: I would like to know exactly, in which aspect blastx and tblasn varies?
As per my understanding, blastx `convert nucleotide to protein and compare with the query protein sequences`
where as, tblastn also does the same thing but, here the query is nucleotide `convert nucleotide to protein and compare with subject protein sequences` Hence, both are doing the same, then what is the difference between these two process?.
Please let me know the difference and why do we use `tblastn` over `blastx` for evolutionary analysis.
Thank you in advance.

Answer: You have mixed-up your query sequences.

Blastx is used when you want to compare a nucleotide sequence to a protein DB. but as we can not compare nucleotides and protein to each other directly the nucleotide sequence is translated (in all 6 frames, this is important !) and then an amino-acid vs aminoacid comparison is done.

For tBlastn, it's the reverse, here you want to compare a protein sequence to a nucleotide DB, and for the same reason as above , the would DB in this case is translated in all 6 frames, to allow protein to protein comparisons.

One could simplify it as: in tblastn, you know what the protein should be like and you want to search for nucleotide sequences where that protein might be similar to. For blastx , you have no idea where or which protein might be present on it so you screen that to a protein DB to look for potential matches of known proteins.

There is no reason why one would prefer one over the other in evolutionary analysis. The only thing is that it's often more convenient to use a protein as query if you want to screen for absence or presence of a certain gene. Moreover evolutionary analyses most often are done on protein level so it's more likely that you want to look for potential other proteins based on a protein you have as query.
Question: I have been creating rainfall plots following the tutorial from [https://bernatgel.github.io/karyoploter_tutorial//Examples/Rainfall/Rainfall.html][1]

Does anyone know how I could extract only mutations with a high degree of spatial clustering from the GRanges object? For example if I want to create a new object which is a list of all mutations that are within 1kb of another mutation. The output would be in this format for input into another program:

Chr start end ref alt
1 10 10 A T
1 50 50 C G
2 181 181 A T
2 280 280 T A

I also have the raw list of all variants in the format above that becomes a GRanges object, so it might be easier just to use a function that works directly on that to find clustered variants? I imagine a tool must already exist to identify variants that are spatially clustered in the genome, but I have not yet found one.

Thanks in advance for suggestions.

[1]: https://bernatgel.github.io/karyoploter_tutorial//Examples/Rainfall/Rainfall.html

Answer: Hi Rubai,

If you are already using `kpPlotRainfall` you can take advatage of the distances it computes. You should call it and assign the result to kp again

kp <- kpPlotRainfall(kp, data=your.data, [...])

and then you can get the distances with

distances <- kp$latest.plot$computed.values$distances

The distances are in log10, so if you want just the variants that are closer than 1kb to the next variant you should

your.data[distances<3,]

And that should do it.

Hope this helps
Question: Hello all,

First post so apologies for any etiquette breaches.

I have a large cohort of over 1000 patients who have had WGS done and then Canvas/Manta structural variant caller applied to their data. I have then merged these into a single VCF file.

In the ID (3rd column) I have the types of call made by one of the two alogrithms e.g. MantaDEL, MantaINV, CanvasGAIN etc etc. After each piece of text there is a string of number representing the graph plot associated with how that call was generated e.g. MantaDel0:1:3:00 which is unique to that call.

I am trying to create individual VCFs of each unique type of call. However, if I run the following code to see what the unique SV types are:

bcftools query -f '%ID/n' SV.vcf.gz

I unsurprisingly get a massive list of each "ID" with the unique number string after it when all I want is the first part, as in to know what are the types of SV call so I can then use bcftools to create individual file types for each.

I hope this makes sense and apologies if this has been covered elsewhere.

Many thanks

Answer: for T in MantaDEL MantaINS MantaBND
do
bcftools view input.vcf | awk -v T=${T} '{split($3,a,/\:/);if(a[1]==T) print $1,$2,$3;}' > ${T}.txt
done
Question: I have a code that is creating me a script where I can change the name of the titles of the sequences in a FASTA file.

This is the text file I'm using:

#Assembly	Genome Center name	RefSeq Accession.version	GenBank Accession.version	NCBI name
GeoFor_1.0	scaffold40	NW_005054297	JH739887	GPS_002009865
GeoFor_1.0	scaffold112	NW_005054298	JH739888	GPS_002009866
GeoFor_1.0	scaffold41	NW_005054299	JH739889	GPS_002009867
GeoFor_1.0	scaffold130	NW_005054300	JH739890	GPS_002009868
GeoFor_1.0	scaffold54	NW_005054301	JH739891	GPS_002009869
GeoFor_1.0	scaffold16	NW_005054302	JH739892	GPS_002009870

This is the FASTA file that I'm using and that I want to change the names. AS you can see, I want to find the scaffold names that match the different JH######.

>Scaffold410	275
TGCATTAATATGAGTGTGTGCTGCAAAAGTTCAGGTCATGGTCCGATCATACTTCACATTTTGGTAGCACTTTAAGCAGAGATCGGTTATCCCATTCTGTGGAAGACTCAACACTATCATAAGGTCCCACAGTTTTATTATCCCTCTGCCTCCCGGAATGCCCCCGGCAGTGAGGGGTACCATCTTCTCAGCAGTAAGGATATTCTTCAGGAGTTCCGTGTGAGCTTTCCCGGATTTAGTTCCATTTTTTAAATACTTCCCAATTCTTTGCTTTG
>Scaffold430	374
CTTTGTTAACTGAAAGAGCCTCTAAGTAGATGACCAGTGCTCAGTTAGTACAGTATGAATTTTGTTTAATGGAACAGGAAGATTTAGTATTGAGAAGCGGTTAAGGGTTTAACCCAGCCTCCTGTCTGAATGGACCTGAAGAGGGGGGCCGGGAAGAAACCCATGACTGCATTAAAGTGATAGATCTCCAGACATGGGCTAGGGAAGATTTACAAGACACTCCCTGGCCTGAGGGAGAAAATATGTTTATTGATGAGTCTTCAAGGGTGGCAGAAGGGAAGCGATTTACAGGATACACAATCATTAATGGAAGGAAATTAAAGGAAGGGGGGAGATTGTCACCCACCTGGTCAGTTCAGACAGCAGAGCTGTAT
>Scaffold1010	597
GGAACACACCTGGGCACACCTGGATGGAGCAGGAACACACCTGGATGGGGTTAGGACACATCTGGATGGCGTTGGGACACACCTGGATGCGCTCAGGGTACACCTG...

Thesis the command I use to create a script to change the names

tail -n +2 scaffold_names_2.txt | while read assemb gcenter refseq genbank ncbi; do echo -ne "sed 's/[[:<:]]$gcenter[[:>:]]/$genbank/g' | " >>script.sh; done

The problem is that I'm not able to save the fasta file with the new names.

This is the last line of my script:

... sed 's/[[:<:]]scaffold4469[[:>:]]/JH767125/g' name.fasta.fa

The script is running without error, but it's doing nothing.

Do you know why? How can I change all the titles and save it as a new fasta file with another name ?

Answer: It's very simple with [SeqKit][1], a cross-platform and ultrafast toolkit for FASTA/Q file manipulation.

Just one command:

cat scaffolds.fasta | seqkit replace --ignore-case --kv-file <(cut -f 2,4 table.tsv) --pattern "^(\w+)" --replacement "{kv}"

Explanation:

SeqKit subcommand `replace` can find substrings with regular expression
and replace them with strings or corresponding values of found substrings provided by the tab-delimited key-value file.

Therefore, the key-value data are firstly extracted by `cut -f 2,4 table.tsv`

$ cut -f 2,4 table.tsv
Genome name
scaffold410 JH739887
scaffold430 JH739888

Then `seqkit replace` command is used to replace FASTA headers. `--pattern "^(\w+)"` refers to the to be replaced substrings (e.g., `Scaffold410`), and the special placeholder symbols `{kv}` represents the corresponding value of the key.

BTW, if you need similar manipulations of CSV/TSV files, [csvtk][2] provides similar functions and features

[1]: http://bioinf.shenwei.me/seqkit/
[2]: http://bioinf.shenwei.me/csvtk
Question: I have assembled bacterial genomes using SPADES. Now I am going to submit them to Genbank, but I need to know the coverage of each assembly. Should I provide the raw read coverage or the filtered final coverage? If the second possibility is true, how do I access these values from the SPADES log file?

Answer: The best way to calculate coverage is by mapping, not by looking at the assembler's logs. For example, with BBMap:

bbmap.sh in=reads.fq ref=contigs.fa covstats=covstats.txt

That will print a message like this:

```
Average coverage: 278.50
Percent scaffolds with any coverage: 100.00
Percent of reference bases covered: 99.98
```

...in addition to creating `covstats.txt` which will list the coverage statistics for each individual scaffold. The reads you use for mapping should be the ones you fed into Spades.
Question: Hello BioStars,

Is there a tool which generates mapping quality score distributions from bam files?

I know I could potentially do this myself, but I am looking for something which would essentially do the work for me, and also something which plugs in nicely to multiqc

I have had a quick search, and it is surprisingly difficult to find a tool which provides this functionality, so I am just wondering if I am missing something obvious here? The obvious candidate would be something like picard, but picard does not seem to offer this functionality from what I can tell

Thanks in advance!

Answer: This may be too simple for what you need, but I have a script [here][1] that uses pysam to extract MapQ scores from a bam file and then generates a kernel density plot of the scores.

Here is the actual code:

#!/usr/bin/env python

import pysam
import seaborn as sns
import argparse
from matplotlib import pyplot as plt

parser = argparse.ArgumentParser(description="Plot mapping quality score from bam file")

parser.add_argument('--bam', required=True, help='path to input bam file', action='store')
parser.add_argument('--png', required=True, help='path to output png file', action='store')

args=parser.parse_args()

output = args.png
input = args.bam

bam = pysam.AlignmentFile(input, "rb")

quals = [read.mapping_quality for read in bam.fetch()]

sns.set_style("dark")
sns.set_context("paper")

plt.figure(figsize=(10,6))

ax = sns.kdeplot(quals, shade = True)
ax.set_ylabel('Density')
ax.set_xlabel('Mapping Quality Score')

plt.savefig(output)

[1]: https://github.com/davidecarlson/Fasta-Manipulation/blob/master/plotMapQ.py