cannin/biostars_qa huggingface dataset:
Question: Hi everybody

I am trying to use htseq for counting the number of reads aligned to genomic features. I use the following command:

htseq-count -s no -t exon -i gene_id --additional-attr=gene_name Alignment-files annotation-files

After running the command the read count table is printed on the terminal screen but there is no spreed sheet saved in my current directory. where can I find the output spreadsheet of htseq-count command?

Best

Answer: Hi Javad,

you need to direct the output to a file :

htseq-count -s no -t exon -i gene_id --additional-attr=gene_name Alignment-files annotation-files > result_file.txt

Cheers,

Michael
Question: Hi, I am running CNVKIT and I am new to it. I was running a batch on two tumor samples with a flat reference.

anaconda3/bin/cnvkit.py batch 1.bam 2.bam -n -t /refseq_exome_genesymbol_headless_CNVkit.csv -f hg19.fasta --access access-5k-mappable.hg19.bed --output-reference my_flat_reference.cnn -d output -p 6

I get

[E::bgzf_read] Read block operation failed with error 4 after 0 of 4 bytes

but it still continues to run. Is it because my bam files are not bzipped? could the results still be usable?

Thanks a lot in advance!
https://www.biostars.org/u/24

Answer: That's odd. CNVkit reads the input BAM files with pysam, which internally uses htslib, the same as samtools. If samtools can read (e.g. index) your BAM files, then they should work fine with CNVkit. I haven't seen that error/warning myself, but if the process continues to run and the outputs from the CNVkit command look reasonable then maybe it's not a problem.
Question: I'm very new to python. I know the basic formula, but I don't really know any commands. I'm learning, but I have to get this done faster than I can learn. I have a tab-delimited file that looks like this:

A	R	R	R	A

R	W	R	W	R

R	M	S	W	K

etc.

and want to generate a new file where I replace certain letters under certain conditions, i.e.:

if all cells in a row are A or R: replace R with G in that row only

if all cells in a row are R and W: replace R with G and W with T in that row only

(+20 other conditions)

so that it looks like:

A	G	G	G	A

G	T	G	T	G

R	M	S	W	K

Nothing should happen if a row has 3 or more unique characters.

The concept would be to read through the first row, see if it matches any conditions, and if it does, write the new output to a file. Iterate through all rows.

file = open(file.txt)
num_columns=numcolumns() # I can manually put in the number of columns for each run, though I'm sure there's a simple command.

for line in file:

if line.contains('A') and line.contains('R') and (count('A')+count('R')==num_columns()):

write.to.newfile.replace('R') with ('G') #I obviously have no idea how to code this

elif:

new_condition #etc. etc. etc.

else:

write line to file as is

Thanks for any help and sorry for the noob question!

Answer: I would probably do it like this. Seems like an odd thing to do, though..care to elaborate about the background?

#!/usr/bin/env python
import sys
import csv
from collections import Counter

# script is called with python scriptname name_of_input_file
# it will create a file named name_of_input_file.out
filename = sys.argv[1]
# with is a safe way to open files for i/o
with open(filename) as fi, open(filename + '.out', 'w') as fo:
# for processing tab-separated, we use csv.reader with delimiter='\t'
for row in csv.reader(fi, delimiter='\t'):
# Counter objects count elements in a collection and store them in a dictionary-like structure
# e.g. {'A': 2, 'R': 3}
symbolCounter = Counter(row)
# this way we can easily determine the number of different characters by checking the length of the Counter object
if len(symbolCounter) == 2:
# assuming all your conditions only take place when there are exactly two different characters
# build a string from the current row
rowString = ''.join(row)
# extract the keys from the Counter object, this will give you the different characters present
# sorting will allow to resolve ambiguities ('AR' == 'RA')
keys = ''.join(sorted(symbolCounter))
# check for your conditions
if keys == 'AR':
# perform the replace-operations on the string built from the row
rowString = rowString.replace('R', 'G')
elif keys == 'RW':
# chaining replace-operations is also possible
rowString = rowString.replace('W', 'T').replace('R', 'G')
# add your other conditions here
# turn the row string back into a list and assign to the original row
row = list(rowString)
pass
# write the row (transformed, if only contains two different characters; untransformed otherwise)
fo.write('\t'.join(row) + '\n')
Question: Hello,

I've got a large number of fastq files generated from a paired end single cell RNAseq experiment.

I'm looking to align them back to mm10 using Hisat2, I can do this if I run every pair individually but is there a way to get hisat to do them all with one code one after the other?

I'm relatively new to the ubuntu coding environment so feel free to dumb it down.
Thanks in advance

Answer: Try this instead and let us know.

for f in `ls -1 *_1.fastq.gz | sed 's/_1.fastq.gz//' `
do
echo hisat2 -x mm10idx -1 ${f}_1.fastq.gz -2 ${f}_2.fastq.gz -S ${f}.bam
done
Question: Does the Common Workflow Language have if statements and for or while loops?
If it does, what is the syntax?

Answer: It does not. This feature is under active discussion at https://groups.google.com/d/topic/common-workflow-language/JU7PSEKr97M/discussion

However, if you have several files that need to be processed (or an array that you want to iterate over for one or more other parameters) then you can use the scatter/gather feature http://www.commonwl.org/v1.0/Workflow.html#WorkflowStep
Question: Does anyone know of a Bioconductor or other R package that can do GO and/or KEGG Gene set enrichment analysis from counts generated with DESeq2? I have a set of counts generated via DESeq2, with some very nice differential expression analysis. I was thinking of redoing the DE analysis using EdgeR which also uses a negative bionomial based normalization, as I could then use the GOana package for GO and KEGG analysis. Can that work with DESeq2, or is there a similar package that can use counts generated in DESeq2?

Answer: GOseq is a suitable R package for GO AND KEGG enrichment. Maybe this link will help you related to DESeq2.

https://wikis.utexas.edu/display/bioiteam/GO+Enrichment+using+goseq
Question: As title.

Because SELEX experiment can detect strong binding site of TF, and ChIP-seq experiment can detect both strong and weak binding site (due to the cross-linked).

So I wonder if PWM derived from SELEX still have the strong binding site property and PWM derived from ChIP-seq still have the strong and weak binding site property.

Or,it will lose the property when these two experiments model to PWM.

Any suggestion or papers will help.

Thanks

Answer: The TFBSs by ChIP-seq data are usually divided into different groups (strong and weak, sometimes with an intermediate group: moderate) based on their binding strength. A common assumption is that "ChIP-seq signal values are correlated with the affinity of the TF-DNA binding" [ [*][1] ]

PWM alone helps us model the sequences at TFBSs but no such information about binding affinity, unless binding scores are also annotated.

[1]: http://www.ploscompbiol.org/article/info%3Adoi%2F10.1371%2Fjournal.pcbi.1003214
Question: Hi,

How to create ribosomal RNA Interval list for to be used with picard CollectRNASeqMetrics?

I have mapped reads to Drosophila_melanogaster Ensembl Genome with Ensembl GTF. I need to identify the percent of reads mapping to ribosomal RNA? Is there a way to create interval list from GTF file?

Thanks!

Answer: The documentation for the IntervalList format is somewhat hard to find. From the linked page:

Represents a list of intervals against a reference sequence that can be written to and read from a file. The file format is relatively simple and reflects the SAM alignment format to a degree. A SAM style header must be present in the file which lists the sequence records against which the intervals are described. After the header the file then contains records one per line in text format with the following values tab-separated: Sequence name, Start position (1-based), End position (1-based, end inclusive), Strand (either + or -), Interval name (an, ideally unique, name for the interval),

So the first thing you need to do is get the header from your SAM/BAM file:

samtools view -H [your.bam] > intervalList.txt

If your GTF file is standard and we assume that it contains only ribosomal intervals, then we need the first, fourth, fifth, seventh, and ninth fields from the file. We can append them onto our text file which contains the header:

cut -s -f 1,4,5,7,9 [your.gtf] >> intervalListBody.txt

This is a very basic approach and you'll probably want to modify it somewhat for your specific needs, but hopefully it's a good start.
Question: Hi,

I would like to filter out some sequences from a fasta file by using a specific pattern.

For example I have this file:

```
>input1
UGAGGUAGUAGG
>input2
CUAUGCUUACC
>out1
UCCCUGAGACCGUGA
>out2
CUCCGGGUACC
>desc1
ACUUCCUUACAUGCCC
```

I know already how I can extract all the fasta sequences with a specific pattern into a new file by using awk.

But what I would like to do is to remove all entries of a specific pattern from the original fasta file and save the newly made file into a new one. In my file above, I would like for example to remove all sequences with the header pattern out. and save only the other to a new file.

Is there a tool somewhere for doing that, or is it possible in awk/sed or even grep

Thanks

Assa

Answer: > I would like for example to remove all sequences with the header pattern out. and save only the other to a new file

awk '/^>/ {P=index($0,"out")==0} {if(P) print} ' in.fasta > out.fasta
Question: I cannot find templates for `phenofile.txt` and `my.phenotype` in the following [PLINK/SEQ code][1]. Can anyone provide the link to the templates? For `phenofile.txt`, is it like the file [here][2]? Can I just have two columns, ID and phenotype? Thanks!

pseq data.vcf.gz v-assoc --phenotype phenofile.txt my.phenotype

[1]: https://atgu.mgh.harvard.edu/plinkseq/assoc.shtml
[2]: https://groups.google.com/forum/#!topic/plinkseq-users/AglI4PWYB00

Answer: The example is like below. You need to load phenotype to your project first.

```
pseq proj load-pheno --file myphe.phe
pseq proj v-assoc --phenotype myphe > output.txt
```

The file `myphe.phe` begins with

```
##myphe
#ID myphe
1 1
2 0
3 1
```
Question: I aligned many fastq files with HISAT2 to grch38. This proceeded without problems.

But in the next step with StringTie, which I am trying to find novel transcripts and their counts with the Gencode27 GTF:

stringtie Donor1_IL2OKT3ZA.HISAT2.sort.bam -G /illumina/runs/RNASeq/Gencode27/gencode.v27.annotation.gtf -A try.tab -p 4 > stringtie.out 2> stringtie.err

However, I get an error

WARNING: no reference transcripts were found for the genomic sequences where reads were mapped!
Please make sure the -G annotation file uses the same naming convention for the genome sequences.

Why doesn't Stringtie recognize Gencode annotation? Do I have to do something to the gencode data?

Update: STAR works with this stringtie, but HISAT2 output doesn't. Strange.

stringtie's output from `cut -f 3 try.tab | sort | uniq`

looks like

703404669@ssxfisctimga004:~/RNASeq_benchmark/GSE96075/HISAT2$ cut -f 3 try.tab | sort | uniq
1
10
11
12
13
14
15
16
17
18
19
2
20
21
22
3
4
5
6
7
8
9
GL000008.2
GL000009.2
GL000194.1
GL000205.2
GL000214.1
GL000218.1
GL000219.1
GL000220.1
GL000221.1
GL000224.1
KI270442.1
KI270706.1
KI270711.1
KI270713.1
KI270721.1
KI270733.1
KI270734.1
KI270742.1
KI270744.1
KI270745.1
MT
Reference
X
Y
chr1
chr10
chr11
chr12
chr13
chr14
chr15
chr16
chr17
chr18
chr19
chr2
chr20
chr21
chr22
chr3
chr4
chr5
chr6
chr7
chr8
chr9
chrM
chrX
chrY

UPDATE: I have found several other instances of this error, but no one ever addressed how to solve this:

https://github.com/gpertea/stringtie/issues/113

https://www.biostars.org/p/247954/

Answer: StringTie gave you an advice there, did you follow up on it?

>Please make sure the -G annotation file uses the same naming convention for the genome sequences.

In case you don't understand that message: it's about the chromosome names. The Gencode annotation uses chromosome names like *chr1*, *chr2* ,*chr3*,... while the grch38 alignments have *1*,*2*,*3*,.. instead. Hence that WARNING message.
Sure for you it might be obvious that "1" is the same as "chr1". But StringTie is made to work on any assembly/genome data, not just the human genome and it's not going to second guess your use of mismatching chromosome names like this..
I would suggest to use the UCSC hg38 genome instead -- or, if you don't want to re-run hisat2 again, find an annotation which uses the 1,2,3.. naming convention for the chromosomes.. or find a *reliable* way to convert the Gencode genomic sequence names to the grch38 naming convention.. This might be as simple as removing the "chr" prefix for the chromosome names -- but this might only work for the main chromosome sequences, the naming convention for the additional/alternate contigs might be different.
Question: Hi,
Is anyone aware of pubilcly available T Cell Receptor repertoire data sets from immuno-sequencing, particularly from tumor or PBMC in cancer patients or mouse. I am looking for a summarized form repertoire profile dataset that includes clonotype (AA and NT), number and frequency of clones, and V, D, J gene name for each clonotype.
Thank you,
- Pankaj

Answer: The paper provided by Denise is about a software for T-cell repertoire post-analysis and thus doesn't contain cancer TCR repertoire datasets (as far as I remember it uses twin repertoire sequencing data from Zvyagin et al PNAS 2014). We have quite a large collection of sequenced repertoires, yet those are mostly from healthy people and people with autoimmune diseases. Regarding cancer datasets you can try to look in these papers:

* http://cancerres.aacrjournals.org/content/early/2012/05/09/0008-5472.CAN-12-0277.full.pdf
* http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4666197/
* http://www.nature.com/articles/srep13664

and check out Immunoseq Analyzer page:

http://ww2.adaptivebiotech.com/IMMUNOSEQ-ANALYZER/DEMO

You'll have to register there and try to find cancer patient samples in published samples section (they can be exported as a tabular file with columns that you need)
Question: Hi,

I'm trying to use a R package for the first time.

Does anyone know how to save/ export a plot (e.g. heatmap) in DiffBind?

Example:
in there example (https://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf)
Heatmap was generated by

plot(tamoxifen)

or

dba.plotHeatmap(tamoxifen)

What should I add in order to save the plot directly from the command? Thanks!!

Answer: The comment above is the correct answer to the question. Before plotting, call one of the functions that redirects the plot to a file, such as `jpeg()`, `pdf()`, or `png()`. After plotting, call `dev.off()` to close the file. You can do multiple plots before calling `dev.off()` and they will all end up in the same file.
Question: Is there a parameter in STAR for setting the minimum length of reads to be mapped?

In the "Log.final.out" file, there is a category, "% of reads unmapped: too short," how do I change the "too short" variable?

Answer: I was able to answer my own question:

The following parameters reduced my "Too short" percent from 10.8% to ~4%. Most of the reads went to multi-mapped or "other," but at least it worked.

--seedSearchStartLmax 30 # increase overall mapping sensitivity
--outFilterScoreMinOverLread 0 --outFilterMatchNminOverLread 0 --outFilterMatchNmin 50 # allow output of short alignments

There apparently isn't a minimum read length threshold for STAR.
Question: Hello, I am analysing ChIP-seq data and I saw the following post by @Devon Ryan (made 1 year, 7 months ago) https://bioinformatics.stackexchange.com/questions/458/when-to-account-for-the-blacklisted-genomic-regions-in-chip-seq-data-analyses/459#459?newreg=dca76bad61c443d7b4f0b1abd1487878 saying that, nowadays with the latest genome assemblies, one has less problems with blacklisted regions since they have been reduced.

I want to know, then, what's the state of the art of this situation? Should I remove them or not? (By the way, I was planning to use Deeptools to do it, but if it's not really necessary anymore I won't). Thank you!

Answer: It's still considered best-practice to remove these regions. For genomes like GRCh38, the blacklisted regions are largely comprised of things like major satellite repeats, which are primarily located in hard-masked telomeric and pericentromeric regions. Given that, these regions will still show aberrantly high signal in all of your samples (thereby skewing normalization and often adding meaningless peaks).
Question: When converting from base count to coverage in [this document][1] they divide by 2.85 billion nucleotides:

```
POP BASE COUNT COVERAGE
ACB 3182037349066 1116.50433300561
```

Does build 37 of the human reference genome still only contain 2.85Gbp? It's the same figure as in [build 35][2]. In my reference sequence I count 3.1Gbp across the autosomes and the sex chromosomes. Have I missed something obvious here?

[1]: ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/sequence_indices/20130502.sequence.index.low_coverage.stats
[2]: http://www.nature.com/nature/journal/v431/n7011/full/nature03001.html

Answer: I have computed this recently on GRCh37 :

(Mappability was computed with GEM program for reads of length 100 and 5 mismatches authorised)

EDIT : And so dividing by 2.85Gb gives you a more realistic estimate of your mean coverage since N's will never be covered by definition.
Question: im using samtool to extract specific regions from the sorted bam file/bam.bai file. e.g if i want to extract the sequence from 329729 to 329840 of chromosome one, then how would i proceed???
the command that im using is

./samtools view alignment.sorted.bam gnl|BGIA2|CA_chr1:329729-329840

and im getting following error:

-bash: BGIA2: command not found
-bash: CA_chr1:329729-329840: command not found
[main_samview] region "gnl" specifies an unknown reference name. Continue anyway.

while the "gnl|BGIA2|CA_chr1" could be found in bam file if is viewed by head command.

further how to extract the multiple sequences with the similar formate in separate file using samtool or any other tool???
thanks

Answer: Hello,

the problem is the `|` in the name. This is interpreted as a pipe to stream the output to the next program, Try to escape this symbol by `\` (Didn't test it)

./samtools view alignment.sorted.bam gnl\|BGIA2\|CA_chr1:329729-329840

Another way could be to create a bed file with the coordinates and use the `-L` parameter.

fin swimmer
Question: I'm using edgeR to test for DE genes in an expression matrix. For each sample if have a condition and a batch. I want to use the glm-functionality in edgeR to test for DE between conditions, while taking into account batches.

Some example data to show my problem. Say if have a count matrix, EM, of 10 sample with the following labels:

```
EM = EM # Imagine matrix of counts here...

conditions = c("con1", "con1", "con1", "con2", "con2", "con2", "con3", "con3", "con3", "con3")
batches = c("batch1", "batch2", "batch3", "batch1", "batch2", "batch3", "batch1", "con2", "con3", "con3")
```

The pipeline could look like this:

```
dge = DGEList(counts=EM) # Create object
dge = calcNormFactors(dge, method='TMM') # Normalize library sizes using TMM
design = model.matrix(~0+conditions+batches) # Create design matrix for glm
colnames(design) = c(levels(conditions), levels(batches)[2:length(levels(batches))]) # Set prettier column names
dge = estimateGLMCommonDisp(dge, design) # Estimate common dispersion
dge = estimateGLMTagwiseDisp(dge, design) # Estimate tagwise dispersion

fit = glmFit(dge,design) # Fit glm

pair_vector = sprintf("%s-%s", "con1", "con3") # Samples to be compared
pair_contrast = makeContrasts(contrasts=pair_vector, levels=design) # Make contrast

lrt = glmLRT(fit, contrast=pair_contrast) # Likelihood ratio test
```

My questions:

1. The design matrix: There is no baseline conditions, so I remove the intersect with the 0+. Is this necessary to do for batches as well, even though they are not directly used as contrasts?
2. Does the glmFit take into account norm-factors for library sizes?

Answer: There are two ways to go about this. Firstly, either allow an intercept and then just use a contrast for the cond3 to cond2 comparison, or don't allow an intercept and group things different (e.g., group=sprint("%s.%s",batch,condition)).

If we just allow an intercept then:

design <- model.matrix(~conditions+batches)
dge = estimateGLMCommonDisp(dge, design)
dge = estimateGLMTagwiseDisp(dge, design)
fit = glmFit(dge,design)
lrt2vs1 <- glmLRT(fit, coef=2)
lrt3vs1 <- glmLRT(fit, coef=3)
lrt3vs2 <- glmLRT(fit,contrast=c(0,-1,1,0,0))

BTW, have a read through the edgeR user's guide, particularly section 3.4.2, which has an almost identical example (it's a really well done user's guide).
Question: Hi,
I'm trying to use sciClone to infer the clonal evolution of some tumor data. I have a couple of questions about its usage (hope I'm not repeating any older question, I checked carefully but none seems satisfying):

1. I have vaf files (from the VarScan variant caller) containing both normal and tumor reads. Does it mean I can create 2 sciClone sample files (1 normal, 1 tumor) from them or should I just use tumor data?

2. What is the meaning of the five columns of the exclude regions file? (I can guess the first three are chromosome, start and end position, not quite sure about the other two)

3. To create the exclude regions file, I can use copy number data from the ICGC releases, where each variant (row) in the file is annotated as "copy neutral LOH", "gain" or "loss". I'm guessing I can populate the exclude regions file with all the rows in the file that are annotated as "copy neutral LOH", is it correct?

I apologize if I'm asking silly questions, I have not much experience in bioinformatics.

Answer: 1) You should be using somatic calls, which means that you should only have one output file. (in your example, created from comparing the tumor and normal using varscan's somatic mode)

2) only three columns are required for the exclude file (chr, start, stop), and that's all that sciClone will look at. One common usage is to call segmented regions of LOH, in which case there may be trailing columns giving info on those segments.

3) if you have copy-number data, you really should be using the segmented copy number data as input if available (which will give the copy number for each segment (which might be not be simple integers of gain or loss, but might be 2.15, or 3.16 etc). If all you have available is the gain or loss designations, then yeah, I guess you could put those into the exclude file. It's not ideal, though.
Question: Hi Biostars,

My question might be a duplicate, but what software would you recommend to simulate RNAseq reads (with different parameters like read length, PE/SE, error-rate, strandedness, etc), by using reference genome and annotations as input.

Google search and omicstools give many results, but I didn't find any soft so far fitting to my needs.
ART looks good (but not specific to RNAseq), and I bet that there would be something in BBtools :)

Would appreciate any comments and suggestions.
Thanks

Answer: `randomreads.sh` in [BBMap suite][1] will do what you want.

I assume you want to generate reads only from the transcriptome (i.e. provide a milti-fasta of transcripts. Don't think of any programs that can take ref genome+GTF file)? Otherwise there is nothing special about RNAseq as far as read generation would be concerned.

[1]: https://sourceforge.net/projects/bbmap/