cannin/biostars_qa huggingface dataset:
Question: SAM/BAM files are that kind of files with hidden treasures that needs to be unrevealed for the inexperienced user. That is me...

My case. After analyzing a BAM file when asking myself how many reads' mates remain unmapped

One possibility to answer this is by analyzing the FLAG values with samtools.

I understand FLAGS are formed by an unique combination of many other individual flags. So, all of these FLAGS values: 73, 89, 121, 153, 185, 137, 77, 141 and so on, contain the "8" , that in turn, should be indicating that the mate read remains unmapped. I got this information from [this WEB page][1] to get an idea about what information the FLAGS can provide

Now a summary..To answer this question I have analyzed a unique BAM file in two ways

- One is by counting the number of `*` present in column 7 (RNEXT value), because in agreement with the official SAM file specification, this could mean that your mate can be unmapped (*This field is set as `*` when the information is unavailable*). In this case, I got over 65000 sequences that could be unmapped
- However, if I run `samtools view file.bam -f 8 | wc -l`, I ended with only 2903 sequences.

One possibility is that when using the `-f` option, the program is looking for a lonely "8" in the FLAG field. But if I look in column 2 in the BAM file, I don't find any FLAG with only that lonely 8. That convinced me that the `samtools view -f FLAG` try to find any combination of FLAG values that intrinsically contains that 8, and thus, it should provide with the information about how many mates are being unmapped

With all this information, I still are not fully confident in knowing what are the right answer to this question. Or I have serious doubts about what is the usefulness of using the `-f` qualificator in the `samtools view` command. Or maybe, many other "lonely" flags should be included in the searching because I miss some important information and/or the orientation do not seem to matter when the BAM file is generated

[1]: http://broadinstitute.github.io/picard/explain-flags.html

Answer: `-f 8` is an `and` operation, so it's just looking that all bits set in `8` are set in the flag (`-F`, on the other hand, just needs a single bit of overlap between the value provided and the flag). Anyway, the most robust method would be `samtools view -cf8 -F 4 file.bam`. This discounts cases where both mates in a pair are unmapped, since I assume you don't care about those. Note that using column 7 to get this is risky since it relies on aligners behaving in a way that isn't guaranteed (e.g., some aligners will give mapping coordinates to unmapped reads).
Question: I am looking for a tool that can randomly shuffle gff features into intergenic regions, but leaving the gene-models 'intact', so that at least all features of a gene (transcript, exon, CDS, etc.) are placed on the same contig and related features are placed *fully inside the interval* of their parent region. Bedtools shuffle doesn't seem to do that, I am trying:

shuffleBed -i genes.gff3 -excl genes.gff3 -g chromsizes.txt -f 0

This command distributes sub-features to different contigs and leads to invalid gene-models, if I add `-chrom`, features are placed on the same contig, but not all features can be placed at all and the resulting gene-models are still not valid. Does anyone maybe have some R-code for this use-case?

Answer: Here's my second attempt, this time in Perl. It does more or less what I want with good speed, maybe I will add a few more options and documentation in the future. This solution is really fast enough, I tried to shuffle 100,000 features thereof 30,000 top-level genes over 30,000 contigs in a few seconds, more or less the time it takes to do the parsing and IO. Feel free to comment on possible improvements bugs, horrific code...

http://mdondrup.github.io/shuffleGff/

Hint, to place genes only in intergenic regions use the same gff file as exclude and input.
Question: I am looking for a database or a way that I could look at broad changes in protein phosphorylation by phosphoproteomics at many targets in a kinase inhibited or a kinase knockdown state. For example, all the proteins that show reduced phosphorylation when Akt (or any given kinase) is inhibited or knocked-down or out.

Any help would be very much appreciated.

Thanks!

Answer: There is no such database for human. Phosphogrid (phosphogrid.org) contains conditional phosphoproteomics data for yeast but even in that case you don't get the actual quantitative scores so there is no easy way to get these types of data right now except to compile them from articles. [Pedro Cutillas][1] has published several studies doing exactly what you want. The two recent papers you might want to have a look are these two:

1. [Large-scale models of signal propagation in human cells derived from discovery phosphoproteomic data][2]
2. [Empirical inference of circuitry and plasticity in a kinase signaling network][3]

There are other studies with fewer inhibitors.We have compiled several of them and we are now trying to publish meta analysis of these types of data. We will make all of the data available under the same protein ID system when that first manuscript gets closer to publication.

[1]: https://scholar.google.co.uk/citations?hl=en&user=bOgIkE4AAAAJ&view_op=list_works&sortby=pubdate
[2]: http://www.nature.com/ncomms/2015/150910/ncomms9033/full/ncomms9033.html
[3]: http://www.pnas.org/content/112/25/7719.short
Question: Here is what I am doing:

#GDS download and processing test
source("http://bioconductor.org/biocLite.R")
biocLite()
library(GEOquery)
gds <- getGEO('GDS756')
#Convert to ExpressionSet
eset <- GDS2eSet(gds,do.log2=TRUE)

eset is summarized like this:

ExpressionSet (storageMode: lockedEnvironment)
assayData: 22283 features, 6 samples
element names: exprs
protocolData: none
phenoData
sampleNames: GSM21712 GSM21713 ... GSM21718 (6 total)
varLabels: sample disease.state description
varMetadata: labelDescription
featureData
featureNames: 1007_s_at 1053_at ... AFFX-TrpnX-M_at (22283 total)
fvarLabels: ID Gene title ... GO:Component ID (21 total)
fvarMetadata: Column labelDescription
experimentData: use 'experimentData(object)'
pubMedIds: 16531451
Annotation:

What I want is to output a table in a text file, containing the following columns: (Affymetrix probe ID, ENTREZ Gene ID, expression for sample 1, ..., expression for sample n).

Weirdly enough, GEOQuery manual says that eset should not contain annotated genes, but it seems from my basic object inspection that the featureData contains just that. However my R knowledge is limited so I was unable to produce my table from it.

Answer: write.table(data.frame(fData(eset),exprs(eset)),file="expression.txt",row.names=FALSE,sep="\t")
Question: Hi,

I'm interested in getting the scientific names of my blast hits ran locally. I see blast+ search apps have option -outfmt which can take sscinames(seems new in version Blast+ 2.2.28), but even using nt from NCBI(no luck with local databases either) I get N\A for this specifier. Similarity for '%S' of -outfmt in blastdbcmd.

For example:

$ blastdbcmd -db nt -entry 229577210 -outfmt '%a || %g || %T || %S || %t'
NM_001743.4 || 229577210 || 9606 || N/A || Homo sapiens calmodulin 2 (phosphorylase kinase, delta) (CALM2), mRNA

Until now I've been using taxids in a very convoluted way. I will get the GIs from my hits, then query the blast db using blastdbcmd to get the taxid and then query the local copy of the NCBI taxonomy database with bioperl to get the scientific name. Now that I see blast+ seems to be able to directly output the scientific name, I would like to simplify things. I'm already able to simplify things a little using the also new output format specifier staxids, so I can now get the taxid directly from the blast output.

So my questions is.

Is there a way to build local blast databases in a way so 'sscinames' can be used to output the scientific name in blast+ results?

In a side note. If there is a way, it seems odd NCBI's nt is not built using it. At least that is the case for the version I got from Jul 11 2013.

Thanks in advance,

Carlos

EDIT: I found I can now use staxids to simplify my life a little. Some additional question formatting. NT updated to version from Jul 11.

Answer: I indeed was able to find my answer in NCBI BLAST documentation:

BLAST Command Line Applications User Manual

Basically the only taxonomy information stored directly in the BLAST database is the taxid. The rest needs to be pulled from an additional database also provided by NCBI:

ftp://ftp.ncbi.nlm.nih.gov/blast/db/taxdb.tar.gz

As per the documentation make sure taxdb database is in the path defined by BLASTDB environment variable. After that you will be able to ask for several additional taxonomic information in the tabular output. For example:

$ blastn -db nt -num_threads 24 -max_target_seqs 1 -outfmt '6 qseqid sseqid evalue bitscore sgi sacc staxids sscinames scomnames stitle' -query 229577210.fasta
gi|229577210|ref|NM_001743.4| gi|229577210|ref|NM_001743.4| 0.0 2418 229577210 NM_001743 9606 Homo sapiens human Homo sapiens calmodulin 2 (phosphorylase kinase, delta) (CALM2), mRNA
Question: I'm interested in getting simple "heterozygous" or "homozygous" designations for all of the samples/SNPs in my multisample VCF file. In the past, I have been using the -GF GT option in GATK's VariantsToTable tool, and then annotating my basecalls in Excel as either heterozygous or homozygous. This takes forever since Excel isn't really built for big data like this. Is there a simple way to output all of the SNPs as 0/1, 0/0, 0/1, or 1/1 instead of C/A, A/A, G/T, C/C? My ideal output would be a txt file in a grid similar to how VariantsToTable outputs data: top row is each sample, while first column is the variant coordinates.

Answer: if you start from a valid vcf file, you can get your desired output simply with this command:

grep -v ^## input.vcf | cut -f1,2,10- | sed 's/:\S*//g'

`grep` to remove all headers but column names, `cut` to select chromosome+position+samples' columns, and `sed` to remove everything but GT from genotype columns.
Question: Dear all,

I am trying to write a fastq file (to be more exact a fastq sanger) from a bam file.
I tried to install bam2fastq (https://github.com/jts/bam2fastq) but I seem to be unable to install it correctly (which bam2fastq on my command line returns nothing). In addition, I tried picard tools, which returns an empty file, even though my bam file is not empty.

1) So I tried to write it via command line, but I have no experience with it, so I am at a loss.

2) In addition to this, I don't understand what the second part of the identifier is. I mean the identifier starts with @ and goes on and then there is a space and there is for example a 1:N:0:ACTTGA. So I wonder, what this 1:N:0:ACTTGA is.

I want something of the form

@HWI-ST999:169:C3001ACXX:2:1101:1193:2063 1:N:0:ACTTGA
NTCCAAAACACTAACCAAGCTTCTTCTTGCTTCTCAAAGCTTTGATGGTTT
+
#1=DDFFFHHHHHJJJJJJJJJJJJIJJJJJJJJJGJJJJJIJJFHIJGHJ

So far I have this:

samtools view input.bam | cut -f1,10 | sed 's/^/@/' | sed 's/\t/\n/' > out.fastq

@HWI-ST999:169:C3001ACXX:2:1101:13677:2145	GATCGGAAGAGCACACGTCTGAACTCCAGTCACACTTGAATCTCGTATGCC

I have been trying this:

samtools view /Users/Mela/Downloads/Galaxy102Unmapped.bam | cut -f1,10 | sed 's/^/@/' | sed 's/\n/\ /+/\ \n/' | head -n 20

Answer: Samtools is not the right tool for this job. Try bedtools bamtofastq (after name-sorting your bam): https://bedtools.readthedocs.io/en/latest/content/tools/bamtofastq.html
Question: Hello,

Thank you for the previous answers to my questions here (https://www.biostars.org/p/135443/#136088). I have a couple of follow up questions.

Once SNPs are discovered, predictions are listed in the fasta files, and they are typically 60bp, with the SNP located in the middle of the prediction. This is true whether I use 100bp reads or 50bp reads. Now, coherence is established by splitting the 60bp prediction into k-mers (k=31) and assesing their k-mer coverage by all samples allowing one mismatch. And if at least one sample has a k-coverage less than 4, the SNP is considered uncoherent.

I was wondering if the length of the prediction can have an effect on coherency? Should I try to lower the k-mer size for shorter reads?

I was also wondering what is the effect of various coverage in different samples? Should we aim to subsample to the same coverage level?

Lastly, I tried genotype discovery with a -P of 6, are the genotypes the ones reported in the fasta file? In some regions I expect more than 2 genotypes.

Thank you,

Adrian

Answer: Hi Adrian.

Without options `-t` or `-T` the length of the prediction is 2k-1 (should be 61 with k=31). With bigger *k*, the prediction is bigger. However, each prediction (without contigs/unitigs extensions) contains exactly *k* k-mers. Thus the number of tested positions for the coherency remains the same with big or small k-mers. By the way, longer kmers may be more difficult to exist in the reads.

I invite you to read the NAR paper in which the effect of kmer lengh is analysed on a set of human reads.

Currently the coverage threshold is applied per read set. However it is not possible yet to fix a coverage per read set. This improvement is necessary and we are currently working on this feature. If your coverages are high and highly dissimilar, it may be a good idea to subsample larger read sets.

Your last question is a bit more tricky. The `-P` behavior depends on the complexity of the combinations of SNPs. Clearly with `-b 0` or `-b 1`, SNPs ---A----C---- + ---A----G---- + ---T----C---- + ---T----G---- for instance are not detected. In this case, with `-b 2` variants are independently detected: (-A----C---/-A----G---), (-T----C---/-T----G---), (---A----C-,---T----C-), and (---A----G-,---T----G-) are detected --> the 2 first detect the C/G SNP with the two possible contexts (A and T) and the two last detects the variant A/T either with context C or G.

Another possibility is ---A----C---- + ---A----G---- + ---T----C---- (without ---A----G-,---T----G-). In this case b 0 detects nothing, while b 1 or b2 detect (-A----C---/-A----G---), and (---A----C-,---T----C-).

I hope this helps.

Best, Pierre
Question: I have a file with the following format, that relates a determined position of the genome to coverage data:

Sequence position coverage

```
NC_006510 1 19
NC_006510 2 22
NC_006510 3 44
NC_006510 4 1
NC_006510 5 1
NC_006510 6 0
NC_006510 7 0
```

... (low coverage positions, under 6)

```
NC_006510 1000 0
NC_006510 1001 66
NC_006510 1002 66
```

I would like to get intervals larger than 800 pb with low coverage along a sequence. For example, a desirable output would be like that:

NC_006510 4..1000 #In this position interval, the sequence has low coverage (<6).

Do you know a clever way to do it?

Thanks in advance.

Answer: If you are looking at 800bp or more stretches with coverage less than or equal to 6X, some thing like following would be useful.

bedtools genomecov -ibam input_sorted.bam -g reference.fa.fai -bg | awk '{ if ($4<=6) print }' | bedtools merge -i - | awk '{ if (($3-$2)>=800 ) print }' | less

Anyway this is not an exact answer to your question, but different approach.
Question: Hi Biostars Users,

I would like to know if you have encountered this issue while working with [CONTRA tool][1] for detecting CNVs in targeted NGS.

[1]: http://sourceforge.net/apps/mediawiki/contra-cnv/index.php?title=CONTRA:_Copy_Number_Analysis_for_Targeted_Resequencing#Examples

In fact, my analysis is "stopping" at the binning process with no error message, I have checked all the required arguments and they seem to be ok :

```
contra apple$ python contra.py --target ngs1_probes.bed --test cic04668.bam --control baseline_ngs1.bed --bed --fasta human_ref.fasta --outFolder ~/contrases

target : ngs1_probes.bed
test : cic04668.bam
control : baseline_ngs1.bed
fasta : human_ref.fasta
outfolder : /Users/apple/contrases
numBin : [20]
minreaddepth : 10
minNBases : 10
sam : False
pval : 0.05
sampleName : No-SampleName
nomultimapped : False
plot : False
bedInput : True
minExon : 2000
largeDeletion : False
Creating Output Folder : Done.
Converting TEST Sample...
DEBUG 123 genomeCoverageBed -ibam cic04668.bam -bga -g /Users/apple/contrases/buf/sample.Genome
Converting CONTROL Sample...
Getting targeted regions DOC...
chr1
chr10
chr11
chr12
chr14
chr15
chr16
chr17
chr19
chr2
chr20
chr3
chr4
chr5
chr6
chr7
chr8
chr9
chrX
Targeted regions pre-processing: Done
Getting targeted regions DOC...
chr1
chr10
chr11
chr12
chr14
chr15
chr16
chr17
chr19
chr2
chr20
chr3
chr4
chr5
chr6
chr7
chr8
chr9
chrX
Targeted regions pre-processing: Done
Test file read depth = 97323185
Control file read depth = 98289477
Pre-processing Completed.
Getting the Log Ratio
Binning....
```

Do you have any idea why?

Thank you for your help

Kiz

Answer: Hi [sruthimouli88][1],

To be honest, I just dropped using CONTRA and I am currently using [ExomeDepth][2] from (R),

I have configured it and it is working really well.. I advise you to use it, and if you need the scripts for CNVs detection, just tell me :)

I hope this helps,

Kizuna

[1]: /u/12128/
[2]: http://cran.r-project.org/web/packages/ExomeDepth/index.html
Question: Hello,

due to "reasons", I have a somewhat important batch effect in a cohort, that is impacting quite negatively my splicing analysis. Is there any way that this could be "dealt with" in KissDE (or with another way ?), like it can be done with DESeq2 (modeling the effect by including it into the design formula : ~ batch + condition) ?

By the way, the new --type1-only parameter in KisSplice 2.6.2 is quite a nice addition, thanks!

Best,

David.

Answer: Hello,

I do not think this is yet possible with kissDE, even though it could be a quick, easy and nice addition to the package! Maybe in the near future...?

Oh, you noticed that option? Cool! Indeed, it can save a lot of time for people only interested in AS events (even though some AS events could be found in the type 4...).

Have you noticed the Shiny Interface that comes with the lastest version of kissDE? We hope that it can help users to quickly navigate and identify interesting entries in kissDE results...

Best,
Audric
Question: Is there any way to convert a cytoband location (e.g., 17p13.3) to its genomic coordinates programmatically (i.e., without having to do a manual conversion)? A web API such as BioMart would be ideal, but I haven't found a way to do it this way. Entrez Eutils would also be okay, but I likewise haven't been able to figure out how to do it this way.

Answer: $ mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -D hg38 -e 'select * from cytoBand where concat(right(chrom,length(chrom)-3),name)="17p13.3"'
+-------+------------+----------+-------+----------+
| chrom | chromStart | chromEnd | name | gieStain |
+-------+------------+----------+-------+----------+
| chr17 | 0 | 3400000 | p13.3 | gneg |
+-------+------------+----------+-------+----------+
Question: I'm using Bowtie2 to align some simulated E coli reads, on my laptop so it's taking a while, and the only measure of progress is the increasing size of the sam file Bowtie is writing. Unfortunately, I don't know how big that sam file should eventually get -- is there any rough guide to estimating the size of a sam file based on the size of the input fastq files? Is it nice and simple and they should be roughly the same size, or will it vary depending on the quality of the alignment?

Answer: Estimation by size is difficult, because reads without mapping will only contribute very little. But you can do something like this:

```
ID=`tail -n1 .sam | cut -f1` # gives you the latest mapped read id
LINE=`grep -n "^@$ID" .fastq | cut -f1 -d:`
# last id line (use > for fasta)
TOTAL=`cat .fastq | wc -l` # total lines
echo "scale=5; $LINE/$TOTAL * 100" | bc # percentage done
```

[edit] the syntax should work now.

[edit] bc scale.
Question: Hello fellow bioinformaticians!

I was wondering if anyone had any clue as to whether or not a tool exists for generating dummy/fake reads from an already assembled genome? So far my searches have yielded very little...

Essentially, the ideal software would take a single-contig assembled genome and generate simulated raw reads based on some input criteria such as desired average coverage, level of base quality, etc... .

People are not always nice enough to release the raw reads from their sequencing projects, hence simulated raw reads is the second best thing.

I realize this might be a very difficult task, but assuming no tool exists for this, do you think it would be possible to make a decent one?

Thank you for the input!

Answer: [BBMap][1]'s 'randomreads' does exactly what you need.

[1]: https://sourceforge.net/projects/bbmap/
Question: I have a table with two columns and I would like to separate its first column in two separated columns or remove some useless values.

Input:

Gene Disease
MSMO1 (6307)	MICROCEPHALY, CONGENITAL CATARACT, AND P... (OMIM:616834)
SC5D (6309)	LATHOSTEROLOSIS (OMIM:607330), Lathosterolosis (ORPHA:46059)
COG1 (9382)	CONGENITAL DISORDER OF GLYCOSYLATION, TY... (OMIM:611209)
COG6 (57511)	CONGENITAL DISORDER OF GLYCOSYLATION, TY... (OMIM:614576)
UBE3A (7337)	15q11q13 microduplication syndrome (ORPHA:238446)

Output:

Gene Disease
MSMO1	MICROCEPHALY, CONGENITAL CATARACT, AND P... (OMIM:616834)
SC5D	LATHOSTEROLOSIS (OMIM:607330), Lathosterolosis (ORPHA:46059)
COG1	CONGENITAL DISORDER OF GLYCOSYLATION, TY... (OMIM:611209)
COG6	CONGENITAL DISORDER OF GLYCOSYLATION, TY... (OMIM:614576)
UBE3A	15q11q13 microduplication syndrome (ORPHA:238446)

Thanks for any solution in advanced.

Answer: library(splitstackshape)
cSplit(df, 'Gene', sep=" ", type.convert=FALSE)
df<-df[-c(column_name_generated_after_split),]