cannin/biostars_qa huggingface dataset:
Question: I started with a 9GB file of a human Xth chromosome with a .bam format, and I am trying to get a specific genic region in .fasta classic format (you know, a one header file, starting with ">", with description on header, and a single line, of continuous tandem single nucleotides).

I have been able to retrieve the reads of the genic region I want in bam format. I share the path I have followed:

```
samtools sort OriginalFile.bam -o OriginalFile.sort.bam
samtools inex OriginalFile.sort.bam
samtools view -b OriginalFile.sort.bam RegName:XX-XX+1 > GeneName.bam
```

after this I got the GeneName.bam file which contains all the reads on which I am interested in, but I can't run the phylogenetic tool I am trying to pipe the GeneName.bam sequence to ... in order to follow my workpath, I need to transform this GeneName.bam file, into this classic .fasta format I described lines ago.

I tried a couple of pipelines in order to achieve this, but what I have got as a retrieved .fasta file, is a man readable file, but with all the single fastq reads, with headers and everything, stacked on top of each other. I need to get the consensus of this, with one header and only contiguous nucleotides.

I got GeneName.fasta file running:

samtools bam2fq GeneName.bam > GeneName.fasta

or

```
samtools bam2fq GeneName.bam > GeneName.fastq
samtools seqtk seq -A GeneName.fastq > GeneName.fasta # (i installed bowtie2, boost, tophat2, seqtk, and bcf tools, as they seem to complement each other's works at some times)
```

and I got the same result on both, analogous file.

With this line of thinking, I thought I may have to run mpileup to the GeneName.bam file, before transforming to fasta, but mpileup gives as a result .bcf format... still, this last assumption is just a form of discussing my problem.

Have anyone out there found the solution for this, can a .bam file converted to classic .fasta format, or can help in any way?

Greetings

Answer: Alright, i solved my problem, so i share what i found to the community so you can solve similar problems

after retrieving a .bam file with the desired length, through the command "samtools view -b ChrName:XX-XX+1 > Outputname.bam"

you need to pipe this output into samtools mpileup command, in order to get a .bcf file

"samtools mpileup -g -f refgenome.fa filename.bam > filename.bcf" (check samtools tutorials to change or check/confirm desired flag options)

then you need to transform .bcf to .vcf callig: "bcftools call -m filename.bcf > filename .vcf"

this last .vcf file needs to be compressed and indexed like following:

bgzip -c filename.vcf > filename.vcf.gz

tabix filename.vcf.gz

after having succesfully run these commands, you can get to the final step, using samtools again:

samtools faidx refgenome.fa ChrName:XX-XX+1 | bcftools consensus filename.vcf.gz > filename.fa

I hope these pipeline is useful for you people, and i thank all the community, and people who aided on the go. Greetings!
Question: Hi all

This question could be considered as a follow up of this discussion.
http://www.biostars.org/post/show/3407/how-to-extract-reads-from-bam-that-overlap-with-specific-regions/#3414
What I need is to extract reads from bam file that fall only within a given region (not overlap the given region), the region being in the form of a gff file or bed file. Overlapping reads could be extracted by several methods (as in the discussion mentioned or BEDTools). The idea is to try to be pretty sure of excluding 5´ UTRs in the process of detecting intergenic transcripts. I saw a tool in BamUtil (http://genome.sph.umich.edu/wiki/BamUtil) called "writeRegion" which would pretty much do what I want. Somehow could not get this running for my dataset.
Was wondering if you guys might have an "R" or some other solution for this.
Thanks in advance
Abi

Answer: I think the options have changed since this reply was written
if you want the default headerless sam file then use the command as per Damian Kao

samtools view input.bam "Chr10:18000-45500" > output.sam

if it's a bam output that you want you will need a '-b'

samtools view -b input.bam "Chr10:18000-45500" > output.bam
Question: Hello,

I'm trying to assemble contigs from a set of paired end reads, denoted as SRR960028_1.fasta and SRR960028_2.fasta. I'm running ABySS on my university's HPC facility. When I run ABySS it terminates early. The line of code itself (within the PBS file) is:

abyss-pe name=abyss_test1 k=63 in='SRR960028_1.fastq SRR960028_2.fastq' v=-v

The tail of the error file looks like this:

Mapped 272979576 of 273907308 reads (99.7%)
Mapped 247491841 of 273907308 reads uniquely (90.4%)
Read 273907308 alignments

Mateless 273907308 100%
Unaligned 0
Singleton 0
FR 0
RF 0
FF 0
Different 0
Total	273907308

abyss-fixmate: error: All reads are mateless. This can happen when first and second read IDs do not match.
error: 'abyss_test1-3.hist': No such file or directory
make: *** [abyss_test1-3.dist] Error 1
make: *** Deleting file `abyss_test1-3.dist'

I've seen the abyss-fixmate error pop up on threads here before, but most of threads about the error seem to have 0 reads in the "Mateless" or the "Total" read section, whereas I have a number. I've also opened the .fasta files and they definitely contain reads. I've seen a few threads that recommend denoting the lines within the .fasta files with a /1 or a /2, but I was under the impression that denoting the files themselves as reads1.fa and reads2.fa would suffice for ABySS (or at least according to the ABySS manual, unless I'm incorrect).

The only thing in the output file is this:

abyss-map -v -j40 -l40 SRR960028_1.fastq SRR960028_2.fastq abyss_test1-3.fa \
|abyss-fixmate -v -l40 -h abyss_test1-3.hist \
|sort -snk3 -k4 \
|DistanceEst -v -j40 -k63 -l40 -s1000 -n10 -o abyss_test1-3.dist abyss_test1-3.hist

The output files generated by ABySS from the run included abyss_test1-1.fa, abyss_test1-2.fa, abyss_test1-3.fa and a abyss_test1-unitigs.fa file. I've checked the head and tail of the files, and they appear to contain contigs.

I'm reluctant to use these files for any analysis because I'm not sure how ABySS assembled them - does anybody know how ABySS assembled them?

Does anybody have any clue as to why ABySS is terminating early, and how I can fix it?

Thanks in advance!

Answer: Hi @jozs2019,

It is most likely a read naming issue that is preventing the reads from being properly paired by ABySS. (I can confirm this if you post the first 10 lines of your read 1 and read 2 files. You can use a gist if you like.)

ABySS requires that the FASTQ IDs (i.e. the first whitespace separated word of the lines beginning with `@`) for reads 1 and read 2 are either identical or have an identical prefix followed by `/1` and `/2`. (See https://github.com/bcgsc/abyss/wiki/ABySS-Users-FAQ#4-my-abyss-assembly-fails-and-i-get-an-error-that-says-abyss-fixmate-error-all-reads-are-mateless-this-can-happen-when-first-and-second-read-ids-do-not-match).

Your results up to `test-3.fa` should be fine, because those first steps of the pipeline don't make any use of the read pairing information. But for the sake of cleanliness and reproducibility, you may want to do a complete rerun of the pipeline after fixing the read IDs.
Question: Hi all,

The documentation for Samtools is minimal at best. I'm still confused on the concept of a clipped read.

1. What is a clipped read? How is it different from a deletion?
2. What is a Soft clip? If the sequence is present in the reference is it different from a mismatch?
3. What is a Hard clip?

Say if I wanted to calculate base pair coverage, would I include soft clipped bases because 'they are present in the ?'

If someone can provide an example such as

```
REF: AGTCG GATCG GTACG
Read: AGTCG xxxCG GTACG
```

That would be even more awesomer

One last question:

I found a ' * ' as a CIGAR string, what does that mean?

Thanks

Answer: Hard masked bases do not appear in the SEQ string, soft masked bases do.

So, if your cigar is: `10H10M10H` then the SEQ will only be 10 bases long.

if your cigar is `10S10M10S` then the SEQ and base-quals will be 30 bases long.

In the case of soft-masking, even though the SEQ is present, it is not used by variant callers and not displayed when you view your data in a viewer. In either case, masked bases should not be used in calculating coverage.

Both of these maskings are different from deletions. Masking simply means the part of the read can not be aligned to the genome (simplified, but a reasonable assumption for most cases, I think). A deletion means that a stretch of genome is not present in the sample and therefore not in the reads.

I'm not sure when H is used instead of the S and vice-versa. I would like to know that.
Question: I'm trying to perform a thorough review of Hi-C and RNAi developments and have found myself looking at publication A (2011) and publication B (2015). B doesn't cite A explicitly but I'm sure there are papers linking the two and I'd like to find them. Are there existing utilities to help with this? I've used NCBI E-utilities in the past but would prefer not to rewrite an application that already exists.

Answer: I played with neo4j and pubmed in 2010: http://plindenbaum.blogspot.fr/2010/02/path-from-egonwillighagen-to-jandot.html but it required to load a bunch of articles.
Question: I just assembled a plasmid using Illumina 2x250 PE. I am pretty confident that the assembly is fine.

I check by mapping the raw data to the assembly. In general I have a very high coverage (sometimes more than 1000x). Also I have very few unmapped pairs.

But I have some regions where the coverage drops to 10-20 fold coverage. This looks concerning, but I still think my assembly is fine because:

1. I this areas I also see barely unmapped pairs

2. I did not see more missmatches from the reads to the assembly sequence as in other regions

Now my questions:

1. What are the properties of sequences where Illumina gives a lower coverage? Any paper? Or a SW to check?

2. How else can I check if my assembly is fine in this low coverage areas?

Answer: Check GC content. Illumina library prep can have problems with both, regions of high as well as low GC
Question: PDB:4DMN is a HIV-1 integrase structure with a bound ligand (PDB:0L9). i want to bring this ligand into alignment with another, reference structure PDB:3AO1.

i can do alignments of the two structures, and i can identify the ligand's heteroatoms by checking the residue.get_full_id() tags. but there are obviously no CA atoms to associated with the ligand to use as references for `Superimposer` to use for alignment. So i figure i need to do it in two steps: first align the two structures, then use the transform found by `Superimposer` on the ligand.

but i don't know how to do the second bit? thanks for any hints,

Rik

Answer: Hi there,

Cross-posting from the Biopython mailing lists.

Superimposer() will give you the rotation/translation matrix you need to superimpose the two structures. Then you just need to apply them selectively to the atoms you want (with Superimposer.apply).

This script posted on [github][1] will do this and an additional sequence alignment step necessary when the two proteins are not exactly the same, in order to get matching atoms to perform the superimposition on.

[1]: https://gist.github.com/JoaoRodrigues/e3a4f2139d10888c679eb1657a4d7080
Question: Hi everyone.

I find some protocol will recommend using

```
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
```

to regress out cell cycle influence or other "noise". But I am confused why we can not drop these genes in `feature selection` steps. For example, we can drop these genes from high variable gene set, so that protecting our latter PCA and cluster analysis from these influence genes.

Best wishes

Guandong Shang

Answer: I still always advise people not to regress out cell cycle status or info. These can be biologically interesting and a big part of many phenotypes. Differing proportions of cycling cells between conditions or samples is both low hanging fruit and something that's easily validated at the bench. I really don't get why people try to remove that, as it's very simple to label populations as "cycling monocytes" and "monocytes" (or whatever) and create supersets as necessary.

I'd prefer that Seurat drop that part of their cell cycle vignette, I have yet to see a case where it's actually helpful.
Question: Hi,

I just downloaded the IGB genomics viewer and would like to use the WS220 version of the C. elegans reference genome. I thought that this version corresponded to the ce10 build according to UCSC:

> This directory contains the Oct. 2010 (WS220/ce10) assembly of the
> C. elegans genome (ce10, Washington University School of Medicine
> GSC and Sanger Institute WS220), as well as repeat annotations and
> GenBank sequences.

Within IGB, when I select the genome version "C_elegans_Oct_2010", the title within the browser window is listed as "C. elegans Oct 2010 (WS140/ce10)".

So now I am confused. Does this genome version within IGB correspond to the WS220 version or the WS140 version of the C. elegans reference genome? Thank you so much for your help!

Answer: Hi smpyonteck,

The C. elegans genome version in IGB "C_elegans_Oct_2010" is the most recent version, i.e. WS220/ce10. The title within the browser "WS140/ce10" is a typo in our synonym file.

Thank you for finding this!

Nowlan
Question: So I have a sam file with extended CIGAR format like this:

17H87=1D12=1D2=3D32=1D4=2D5=2D13=....

And I want to convert it to a sam file containing regular CIGAR format that will look like this:

17S83M1D14M1D4M3D31M1D2M2D6M...

(they do not represent the same reads... I want only to visualize the problem)

Anyone knows some software/script that will convert sam/bam with extended CIGAR to a sam/bam with regular CIGAR?
It seems like a standard problem that might already have a solution so I want to ask here before solving it the hard way(coding)

Answer: I've quickly written one: https://github.com/lindenb/jvarkit/wiki/Biostar234081

$ cat toy.sam
@SQ	SN:ref	LN:45
@SQ	SN:ref2	LN:40
r001	163	ref	7	30	1M2X5=4I4M1D3M	=	37	39	TTAGATAAAGAGGATACTG*XX:B:S,12561,2,20,112

$ java -jar dist/biostar234081.jar toy.sam
@HD	VN:1.5	SO:unsorted
@SQ	SN:ref	LN:45
@SQ	SN:ref2	LN:40
r001	163	ref	7	30	8M4I4M1D3M	=	37	39	TTAGATAAAGAGGATACTG*XX:B:S,12561,2,20,112
Question: Dear all,

I am testing TAXAasign and getting the following error:

[esther@localhost testes]$ ~/Softwares/TAXAassign-master/TAXAassign.sh -t 70 -m 60 -a "60,70,80,95,95,97" -f test.fasta
[2017-05-02 17:10:09] TAXAassign v0.4. Copyright (c) 2013 Computational Microbial Genomics Group, University of Glasgow, UK
[2017-05-02 17:10:09] Using ~/ncbi-blast-2.6.0+-src/c++/ReleaseMT/bin/blastn
[2017-05-02 17:10:09] Using ~/Softwares/TAXAassign-master/scripts/blast_concat_taxon.py
[2017-05-02 17:10:09] Using ~/Softwares/TAXAassign-master/scripts/blast_gen_assignments.pl
[2017-05-02 17:10:09] Blast against NCBI's nt database with minimum percent ident of 60%, maximum of 10 reference sequences, and evalue of 0.0001 in blastn.

BLAST Database error: No alias or index file found for nucleotide database [~ncbi-blast-2.6.0+-src/c++/ReleaseMT/bin/nt] in search path [~/testes::]

I followed all steps here (https://github.com/umerijaz/TAXAassign). When I run update_blastdb.pl nt the tar.gz files are saved at: ncbi-blast-2.6.0+-src/c++/ReleaseMT/bin
I created the folder nt inside ncbi/bin and moved the tar.gz files there. I also tried to untar. Not working yet.

Any idea?

Kind regards

Answer: Have you tried to set BLASTDB variable to the location where the nt files are located?
Question: Hello,

I have aligned the RNA-seq data to the human genome and I used FASTQC for the pre-processing data. Can I use FASTQC for post alignment QC also? Or is there a better way of doing it?

Thank you in advance.

Answer: RSeQC is used a lot, but I find [QoRTs](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4506620/) a bit more user-friendly, especially if you have numerous samples.

In addition, feeding the results of a) FastQC, b) STAR (or whatever aligner you've used), and c) featureCounts to [MultiQC](http://multiqc.info/docs/) is already quite useful.

The typical things you want to look out for:

* at least 80% alignment rate
* not too many intronic/intergenic reads
* even gene body coverage
Question: Hi to all.

I have a really basic question. Studying ChIP coupled with q-PCR or sequencing, I see different levels of a certain protein among a genomic region. I cannot really understand how can we have different protein levels. I thought that the binding of a protein in a region is bivalent-it occurs or not. But, as I can understand, we can have different concentration/quantity of protein bound to a region? And if this is the case, how in ChIP-qPCR we get finally DNA fragments, so how can we measure this binding from the quantity of fragments?

Thank you very much

Answer: While a protein is either bound or not to a particular chromosome in a particular cell, you're not doing ChIP qPCR on single cells, but rather on a bulk of cells. Heck, even within a single cell you can have different amounts of binding due to the cell's ploidy.
Question: Hi,

I just started to work with single end reads, which are already trimmed for adapter sequences and quality.
Do I have to trimm the reads now to the same length of e.g. 100nt for mapping them with STAR?
Is there a negative effect, if I don't?

Answer: If they are already trimmed for adapters and quality, don't trim more. Trimming will make sequences shorter, and shorter sequences tend to map more to multiple locations.

What is the length range of your reads? I generally keep reads only within a certain range, and discard the shorter reads. For example, for a 100bp dataset, I keep reads from 70-100bp after trimming, and discard the rest.
Question: I have a salmonella treated vs KD condition , i want something like this but I not sure how linear regression can be done

Gene	SCRSAL1	KDSAL3
CXCL1	10.8273907686	11.2751328718
CXCL2	8.001064472	8.5753884341
CXCL3	8.1965832765	8.7086794002
CXCL5	8.0973357585	8.6242736829
CXCL6	9.9053940183	10.8940613053
CXCL8	13.1127454083	13.4402538026
CXCL10	7.2590758038	7.8202130661
CCL20	9.0933584064	9.9218889828
CCL22	7.3643976273	7.7911788897
CCL28	8.3723507486	8.6786344048
IL36G	4.3423919715	4.5673548181

This is my dataset the kind of figure i m trying is this ![enter image description here][1]
Any suggestion or help would be highly appreciated

[1]: https://i.imgur.com/72t9hkq.png

Answer: With perfect data, if you plot gene levels for knockdown vs control, any gene above the 45˚ line is upregulated in the knockdown compared to control, the farther away from the line, the stronger the effect. A linear regression will help in dealing with imperfect data. To help decide on which genes are upregulated, you can calculate (and plot) a confidence interval on each side of the regression line.
For how to do this in R, see this [Rstudio notebook][1].

[1]: https://rpubs.com/aaronsc32/regression-confidence-prediction-intervals
Question: I have several datasets in vcf format(about 1GB per file) and I want to extract snp info in each samples. I usually use python so I tried pyvcf, but it works slow in such large datasets. Are there any better ways to do that, especially in python? Thanks in advance

Answer: A trick is to compress the VCF file with bgzip, and index it with tabix. Both these tools are downloadable from the tabix home page.

See also this page for more documentation. Follow the instructions there:

bgzip my_file.vcf
tabix -p vcf my_file.vcf.gz

I am not sure if pyvcf supports compressed and indexed files. According to the source code, it seems so.

EDIT: yes, it seems that pyvcf supports compressed and indexed VCF files. See the example in pyvcf documentation:

>>> vcf_reader = vcf.Reader(filename='vcf/test/tb.vcf.gz')
>>> # fetch all records on chromosome 20 from base 1110696 through 1230237
>>> for record in vcf_reader.fetch('20', 1110695, 1230237):
... print record
Record(CHROM=20, POS=1110696, REF=A, ALT=[G, T])
Record(CHROM=20, POS=1230237, REF=T, ALT=[None])
Question: Hello folks!

I want to mask a genome with a particular repeat library using RepeatMasker.

Then I want to cross the coordinates of the repeats with those of gene annotations to find overlaps between them and study associations and stuff.

I'm only starting to consider feasible ways to do that so any input would be great.

Thanks!

Answer: Sounds like a job for [BedTools][1]. You should be able to make a GFF from the RepeatMasker output and use `bedtools intersect` to find overlaps.

[1]: http://bedtools.readthedocs.org/en/latest/
Question: Dear all,

I have pair-end RNA-seq data (Illumina) from parasite and I would like to do De-Novo assembly by TRINITY. I have reference genome of my host organism so I can map my data to host and remove from fastq contaminations.

My plan is:

1. Map with bwa/bowtie/novoaling my pair-end FASTQ files to a host reference genome
2. Remove hits from fastq files (cleaning contaminations)
3. For the rest of FASTQ files use TRINITY for De-Novo transcript assembly

My question is:

May I use aligners (bwa etc.) and align raw fastq files to host DNA and then remove contaminants from fastq files? Question is because my data are from RNA-seq project NOT DNA.

How can I remove the sequences from raw fastq files that align to host DNA (cleaning process)?

Or if you have any other advice how to prepare data to TRINITY pipeline I will appreciate it.

Thank you so much for any comment and sharing your experience.

Answer: If you have RNAseq data, you'd be better to stick with an aligner intended for spliced alignments (e.g. STAR). Most of these have an option to place unmapped reads/pairs in a new fastq file(s), which you could then feed to trinity or any other assembler (i.e, step #2 will be done for you). I don't have any advice on good assemblers, hopefully others will chime in with feedback there.
Question: I am attempting to call mutations of raw RNA-Seq of single prostate circulating tumor cells (paper) without matching normals.

From cursory research I found that tumor only without a matched normal calls 117x more variants than matched tumor/normal pairs (Appistry). I also found an interesting study that (successfully) discriminated somatic and germline mutations without matching normals but using virtual normal (i.e., normal samples from unrelated individuals) (paper).

My question is, are blood tumor samples treated any differently than tissue tumor samples in respect to calling variants? If there are no normal blood samples available, is filtering against the EVS dataset my best option? Thank you for your time and help.

Answer: Some of this has been covered in previous questions on here (dig around a little), but the gist is this

- you will never get rid of all germline mutations without a matched normal.
- make sure your variant caller is appropriate and doesn't bias you towards mutations at 50%/100% VAF, as some germline callers can
- a panel of normals can be very helpful in filtering out both common population variants and sequencing artifacts. Every individual contains private rare mutations that will ultimately be indistinguishable from a somatic hit, though.
- yes, use EVS (or maybe even the non-TCGA-derived part of ExAC) to remove common variants.
- RNAseq complicates things even further, because you're dealing with an enhanced error rate.