cannin/biostars_qa huggingface dataset:
Question: Hi !

I started using snakemake to replace my bash scripts in order to have a more "pretty" code, but I have some problems, especially to parallel my jobs.

I have this code:

FILES = [ os.path.basename(x) for x in glob.glob("Experience/*") ]

SAMPLES = list(set([ "_".join(x.split("_")[:2]) for x in FILES]))

CONDITIONS = list(set(x.split("_")[0] for x in SAMPLES))

for path in DIRS:
if not os.path.exists(path):
os.mkdir(path)

rule all:
input:
expand('Trimming/{sample}_R1.trim.fastq', sample=SAMPLES)

rule trimming:
input:
adapters = ADAPTERS,
r1 = 'Experience/{sample}_R1.fastq.gz',
r2 = 'Experience/{sample}_R2.fastq.gz'

output:
r1 = 'Trimming/{sample}_R1.trim.fastq',
r2 = 'Trimming/{sample}_R2.trim.fastq'

message: ''' --- Trimming --- '''

shell: ' bbduk.sh in1="{input.r1}" in2="{input.r2}" out1="{output.r1}" out2="{output.r2}" \
ref="{input.adapters}" minlen='+str(minlen)+' ktrim='+ktrim+' k='+str(k)+' qtrim='+qtrim+' trimq='+str(trimq)+' hdist='+str(hdist)+' tpe tbo '

I have 5 samples, having used the wildcard "sample", I was expecting that my 5 trimming start at the same time, but they start one after the other ..
What's wrong with my code?

thank you in advance

Answer: How did you start the pipeline? Did you use the `-j` option of snakemake to allow multiple jobs?
Question: Kind of a basic awk question. I have a bed file where I have been doing the following to limit interval sizes to those under 1000bp with the following awk script, how would I go about doing the same thing but limiting the output to those between 150-200bp only? Thanks!

awk '{if($3-$2 <= 1000) print}' test.bedpe > test_under1000.bedpe

Answer: Use awk `AND` operator:

awk '{if (CONDITION && CONDITION) print}' test.bedpe > test_between_150-200.bedpe
Question: Hi everyone,

I am trying to make a consensus from a aligned bam file.
This is my workflow:

bcftools mpileup -Ou -f reference_HG19.fa aligned.bam | bcftools call -mv -Ob -o calls.bcf

bcftools index calls.bcf

cat reference_HG19.fa | bcftools consensus calls.bcf > consensus.fa

However, the created consensus.fa file turned out to be identical with the reference_HG19.fa file.
I looked at calls.bcf file with bcftools view, and it seems the following columns are all empty:

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT aligned.bam

When I look at calls.bcf with samtools view I get the error "Aborted".

Does anyone have an idea how to solve this issues?
Thank you very much!

Robert

Answer: We finally found the solution:
The downloaded `bam` was generated by paired end sequencing. For some unknown reason `bcftools` considers these read pairs to be "anomalous". By using the `-A` input option one prevents `bcftools` from skipping anomalous read pairs. Now it's working as it is supposed to.
Question: Dear all,

referring to the batch correction methods for scRNA-seq, would you have any preference and/or comments ? among possible choices :

-- MNNCorrect, as outlined in SimpleSingleCell workflows :

https://bioconductor.org/packages/release/workflows/html/simpleSingleCell.html

-- ZINB-WAVE :

https://bioconductor.org/packages/release/bioc/html/zinbwave.html

-- HARMONY :

https://www.biorxiv.org/content/10.1101/461954v2

-- SCTransform :

https://satijalab.org/seurat/v3.0/integration.html

thanks a lot,

bogdan

Answer: So far, MNN is the best (but still very limit) algorithm for general batch effect correction method. But based on the recent paper (https://www.nature.com/articles/s41587-019-0113-3 ), in some situations, it just exhibits minor improvement than doing nothing. It's all depends on how good your data are.
Question: ***EDIT SKB 17MAR16***
Hi Group,

figured out that I was not calling the correct variable in my loop. I made a handle file outside of the loop to append, but then was using the SeqIO.write function without the handle, so it was writing over the same file with just the last record. Please find the correct script below and thanks for the help!

# this script is used to convert fastq files to fasta files
# then to rename the fasta ID with the sample ID from the lab

from Bio import SeqIO
import sys

# grabbing the file and the name
seq_file = sys.argv[1]
labels = seq_file.split(".")

# converting the file from fastq to fasta
SeqIO.convert(seq_file,"fastq",labels[0]+".fasta","fasta")

# taking the converted file and then changing the fasta header
handle = open(labels[0]+".fasta","a")

for seq_record in SeqIO.parse(handle,"fasta"):
old_header = seq_record.id
new_header = labels[0]
seq_record.id = new_header + "_" + old_header # renaming the pseudogene with
# the lab id and the referance
# used
seq_record.description = "" # this strips the old header out
SeqIO.write(seq_record, handle,"fasta")

handle.close()

***/EDIT***

Hi group,

I wrote a script that will take the fastq output from samtools mpileup > vcf2fq and create a fasta file that includes the drops in coverage. When I use this scripted for a bacterial genome it works as expected, however when I use it on a segmented genome it only outputs the last segment in the file. I'm using python3.4 and biopython1.66. Please see my code and troubleshooting below.

# this script is used to convert fastq files to fasta files
# then to rename the fasta ID with the sample ID from the lab

from Bio import SeqIO
import sys

# grabbing the file and the name
seq_file = sys.argv[1]
labels = seq_file.split(".")

# converting the file from fastq to fasta
SeqIO.convert(seq_file,"fastq",labels[0]+".fasta","fasta")

# taking the converted file and then changing the fasta header
handle = open(labels[0]+".fasta","a+")

for seq_record in SeqIO.parse(handle,"fasta"):
old_header = seq_record.id
new_header = labels[0]
seq_record.id = new_header + "_" + old_header # renaming the pseudogene with
# the lab id and the referance
# used
seq_record.description = "" # this strips the old header out
SeqIO.write(seq_record, labels[0]+".fasta","fasta")

handle.close()

I think I have an error when writing the new record. I have been playing with this line:

handle = open(labels[0]+".fasta","a+")

When I have it written like this, It runs without error, I get a multi-fasta output without the header change; however when I have the code written like this:

handle = open(labels[0]+".fasta","rU")

like in the biopython tutorial, I get the output header I want, but only the last record in the file. I'm not sure if or why with the file open in "append" is just skipping the next part of the script. I would appreciate any input be appreciated

Answer: See original thread for answer
Question: Hello, I am trying to extract some transcript sequences from a stringtie merged gtf using gffread and am getting the following errors:

Error (GFaSeqGet): subsequence cannot be larger than 227
Error getting subseq for MSTRG.17.1 (1..229)!

This error happens for many of gene entries, some but not all transcripts end up getting extracted. The gtf file was created in StringTie using the same reference genome file as everything else (where my bam files are from, the reference annotation gtf file). This annotation is from running StringTie first and merging the gtf annotation files. The main problem seems to be that the end coordinates exceed sequence length when extracting transcripts even though I am sure that they were made correctly. If there is a way to just skip that particular gene that could work as well.

This is the command I used:
gffread /Users/Desktop/StringTie_comp_annotation_new.gtf -g /Users/Desktop/Xla.v91.repeatMasked.fa -w transcripts.fasta

Answer: Hi, jjrin.

I'm having a similar situation here. I have a BAM file from Bowtie2 and a BAM file from BWA. The Bowtie2 BAM file works fine but when I run gffread with the BWA BAM file I got the same error. I found in other threads that it may be related with soft-clipped reads extending beyond the contig.

Anyway, I got a workarround by just replacing the end position of the GTF file by the sequence size on my reference FASTA file. If this solution is a valid one to you, I let the script here: [gtf_fixer_togffread.py][1]

[1]: https://github.com/VictorGambarini/gtf_fixer_to_gffread
Question: Hi All

I have some samples with low mapping in Salmon (40% and less) that have higher alignments in Tophat, and trying to troubleshoot.

I picked some of the unmapped reads (from writeunmapped salmon parameter) and Blat them to human.

Some have 2 or more matches with identity 99% to 100%
And some have many many matches, I need to scroll the page down too much. Many of these matches are 100% and some range between 85% to 100% identity.

I looked also into the “ambig_info.tsv’ , found some records with 0 unique mapping and more than 100 ambiguous mapping, but couldn’t relate to those unmapped.

This is how one match of one of them look in one mate and the other:

ACTIONS QUERY SCORE START END QSIZE IDENTITY CHRO STRAND START END SPAN

browser details YourSeq 22 100 122 151 100.0% 10 + 37378275 37378303 29
browser details YourSeq 22 52 74 151 100.0% 10 - 37378275 37378303 29

So why this is not counted as mapped for example? Any hint, clue?

Thanks

Answer: Salmon maps reads to the transcriptome while TopHat aligns to the genome (assuming that you are using Tophat in the "normal" way). I suspect that if you examine the reads that map with TopHat but not with Salmon, some (perhaps large) proportion of them is mapping to intergenic and/or intronic regions, so those reads will not be mapped by Salmon. Taking a look in a browser at the BAM file from TopHat or using a tool like RNASeqQC can help quantify the intergenic/intronic reads.
Question: Hello world !

I am working on microarray data (Hugene 2.0 st from affymetrix) and I was willing to get probe-level informations.

For one probeset (fsetid), I need all the different probes (fid) and their chromosome, start, stop, sequence etc.

I know I can get this on NetAffx but as I have to check multiple probesets, I wanted to use the library *oligo* in association with *pd.hugene.2.0.st*

Here is my code tested on a 'test-list' of 6 .CEL files:

## Library loading
library(oligo)
library(dplyr)
library(pd.hugene.2.0.st)

# .CEL importation and read
celFiles <- list.celfiles('./', full.names=TRUE)
batch = read.celfiles(celFiles)

# Expression matrix tmp
tmp=batch@assayData$exprs
dim(tmp)

# Connexion pd.hugene.2.0.st@getdb() and normalization
conn = pd.hugene.2.0.st@getdb()
norm = rma(batch, background=TRUE, normalize=TRUE, target="core")

# Data importation
tableNames=dbListTables(conn)
tablesList=list()

for (tableName in tableNames){
tablesList[[tableName]]=dbGetQuery(conn,paste('SELECT * FROM',tableName))
print(tableName)
str(tablesList[[tableName]])
}

# Simplification of the tablesList dataframe
testList <- tablesList$featureSet
testList2 <- tablesList$pmfeature
testListComplet <- dplyr::full_join(testList, testList2, by = "fsetid")

# Test on one probeset randomly selected : probeset 16730541
ind = testListComplet$fid[which(testListComplet$fsetid == 16730541)]
tmp[ind,]
probeset_test <- subset(testListComplet, fsetid == 16730541)
probeset_test

Here is the probeset_test result I get :

fid fsetid strand start stop transcript_cluster_id exon_id crosshyb_type level chrom type fid atom x y
153940 16730541 0 102218019 102218094 16730540 5033332 1 NA 11 1 944681 153940 48 586
153941 16730541 0 102218019 102218094 16730540 5033332 1 NA 11 1 2165196 153941 279 1343

So, according to pd.hugene.2.0.st, for probeset **16730541**, I have 2 probes : **153940 & 153941** with a common start (102218019) and stop (102218094). But when I check on the NetAffx website ( https://www.affymetrix.com/analysis/netaffx/exon/wtgene_probe_set.affx?pk=712:16730541 ), here are the two probes I get :

atcagcggcgccgacaaggagatac	chr11:102218019-102218043 (+)

cagcaaacacggaagctgcgcggct	chr11:102218070-102218094 (+)

Did I do something wrong ? Did I misunderstand something ?
Could you please enlighten me ?!

I am sorry if this is a stupid question but I could not figure it out by myself :(

Thank you so much in advance.

Here is my R info :

platform x86_64-apple-darwin15.6.0
arch x86_64
os darwin15.6.0
system x86_64, darwin15.6.0
status
major 3
minor 4.0
year 2017
month 04
day 21
svn rev 72570
language R
version.string R version 3.4.0 (2017-04-21)
nickname You Stupid Darkness

Answer: You are right (almost!). If you see further in the NetAffx page that you linked

**Probe Set Location:** The genomic location of the probe set for the genome assembly used at array design time (See Genome Source). These coordinates ***begin at the first base of the first probe sequence and end at the last probe of the probe set.***
Question: Hi

I know by

sp_url <- paste("https://cancer.sanger.ac.uk/cancergenome/assets/",
+ "signatures_probabilities.txt", sep = "")

implemented in `MutationalPattern` R package

I can download mutational signatures from the COSMIC website but version 2 (30 signatures)

Anyone knows where I can obtain these probabilities for version 3 (60 signatures)

Thanks in advance

Answer: [Here][1] you can download SBS, DBS and ID probabilities for PCAWG reference signatures. Make sure which one are you using, there is one from SigProfiler and the other from SignatureAnalyzer, two methods used in original publication.

[1]: https://www.synapse.org/#!Synapse:syn12009743
Question: Hi Guys,

I have a BAM file, and a big read list. What I want to do is to remove the reads in the read list from the BAM file. I can transform Bam to Sam file and then use a Python script to remove unwanted reads. And then transform Sam to Bam again. But I am wondering if there is a more efficient way, which I mean faster, easier, and memory-efficient, to achieve this goal?

Any advice is appreciated!

Tao

Answer: [picard FilterSamReads](http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads)

> **READ_LIST_FILE (File)** Read List File containing reads that will be included or excluded from the OUTPUT SAM or BAM file. Default value: null.
Question: Hi.

I post a question continuously.

In this question , curious things are about somatic and germline.

Please see the below paragraph with bold text.

At every position where one or both samples had a variant, VarScan performs a direct comparison between tumor and normal genotype and supporting read counts to determine the somatic status. (I do understand so far but next sentence is confused to me)

**Variants present in both samples are classified as somatic , variants heterozygous in the normal but homozygous in the tumor are classified as LOH(loss of heterozygosity) , and variant shared between samples are classified as germline**.

the second sentence above, What is the difference between somatic and germline?

The sentence that Variants present in both samples are confused in meaning. I knew that somatic mutation is that the variants only exist in tumor sample, but reading and thinking again that sentence .., I can't connect the somatic concept I already knew.

and also I post a question before , but re-question about homozygous and heterozygous.

how can we define homo or hetero considering we have a only one pair of chromosome ?

I can't imagine the homo or hetero in next-seq, because I knew that the only one pair of chromosomes is sequenced...

Am I misunderstanding about concept?

I hope you understand my fuzzy questions.

I am looking forward to your answer

Thank you!

Answer: A somatic variant (=mutation) is a mutation that happens at some point in your life. This variant differs to what is hardcoded in your genome at birth. In cancer genomics, e.g. if you want to find somatic mutations, which means acquired mutations that are unique to the cancer cells (either as a cause or a consequence of the cells transforming into tumor cells), you usually compare the tumor genome against the wildtype (germline) genome of the same individual (e.g. PBMCs, healty tissue). If you variant in the tumor but not in the germline, this could be a hint for a somatic mutation.

Here is a simplified example explaining hetero and homozygosity:

Lets assume, we have a diploid cell (which is not always the case for e.g. tumor cells which can have copy number variations). If variant A (e.g. a SNV) is homozygous, you expect to see 100% of the reads from your exome-seq or wgs (lets assume NGS doesnt produce errors) containing this variant. If this variant is heterozygous (so one chromosome has the variant A and the other has the wildtype B), you expect to see 50% of the reads to contain A and 50% of the reads to contain B.

Note, that this explanation is for DNA-Seq, not for RNA-Seq.

Maybe this publication is also useful for understanding the very interesting field of allele frequencies:

http://www.nature.com/srep/2014/140422/srep04743/full/srep04743.html
Question: How can I create tabix .tbi index file for some .vcf file using java (htsjdk)?

Answer: Use a [TabixIndexCreator][1]

[1]: http://samtools.github.io/htsjdk/javadoc/htsjdk/index.html?htsjdk/tribble/index/tabix/TabixIndexCreator.html
Question: Hello,

I would like to convert a blast tabular format to bed format, I am interested in the subject ID, start and end, so columns 2, 9, 10. Hence it seemed to me this simple awk would do:

awk '{print($2"\t"$9-1"\t"$10)}' file.blastn > file.bed

However, when attempting a bedtools intersect

bedtools intersect -a file.annotation.gff -b file.bed
Error: unable to open file or unable to determine types for file file.bed

- Please ensure that your file is TAB delimited (e.g., cat -t FILE).
- Also ensure that your file has integer chromosome coordinates in the
expected columns (e.g., cols 2 and 3 for BED).

I did check that the file is indeed TAB delimited.
However, I realised that the coordonates difference between columns 10 and 9 was sometimes negative. Example

awk '{print $3 -$2}' AZF180.bed|head
-200
-214
-253

It's not clear how I should deal with that. Simply swapping the columns such as rows with negative sum values become positive?

Thanks for the insight

Answer: I'll assume you have tabular blast output (if not then you should re-run and ask for tab output of blast)

In the tabular output if the alignment is forward against reverse the coordinates will indeed be large to small (instead of small to large). To have correct bed format you will need to switch those coordinates around indeed.

this can be achieved in the same awk cmdline you already used, expanding it a bit (tip: if-condition)
Question: I'm having trouble with Bedtools Intersect. It is not reporting overlap when all of feature A is within feature B.

Command:

bedtools intersect -a File1.bed -b File2.gff -wao -f 0.90 -sorted > Out.file

Out.file:

VMNF01000005.1	484754	484980	.	.	+	.	.	.	-1	-1	.	.	.	.	0
VMNF01000005.1	484754	484980	.	.	+	.	.	.	-1	-1	.	.	.	.	0
VMNF01000005.1	484754	484980	.	.	+	.	.	.	-1	-1	.	.	.	.	0
VMNF01000005.1	484754	484980	.	.	-	.	.	.	-1	-1	.	.	.	.	0
VMNF01000007.1	6294425	6294650	.	.	-	.	.	.	-1	-1	.	.	.	.	0
VMNF01000007.1	6294425	6294650	.	.	-	.	.	.	-1	-1	.	.	.	.	0
VMNF01000007.1	6294425	6294650	.	.	-	.	.	.	-1	-1	.	.	.	.	0
VMNF01000007.1	6294425	6294650	.	.	+	.	.	.	-1	-1	.	.	.	.	0
VMNF01000008.1	1441418	1441616	.	.	-	.	.	.	-1	-1	.	.	.	.	0

When I search manually, all of feature A sits within feature B.

![Image of overlap when viewed in Artemis. Feature A (pink and yellow highlighted region within circle) sits within feature B (red highlighted region).][1]

Can anybody please explain where I have gone wrong?

Bedtools version: bedtools v2.29.2

[1]: https://i.ibb.co/gy5x7V5/Artemis-VMNF-Promoter-in-red-mimp-in-pink.png

Answer: I will leave this for those who may be interested - the issue was that the scaffolds/contigs had a slight difference in the names between the .bed file and the .gff file. The .bed file included the suffix ".1". After adding ".1" to the scaffold/contig names in File2.gff, both bedtools intersect and bedmap worked.

E.g.

=> File1.bed <==
VMNF01000005.1 482754 486980 . . +
VMNF01000005.1 482754 486980 . . +

==> File2.gff <==
VMNF01000002 GenBank PROMOTER 1 34 . + 1 ID=FocTR4_00017231;Name=FocTR4_00017231
VMNF01000002 GenBank PROMOTER 4227 5039 . - 1 ID=FocTR4_00017232;Name=FocTR4_00017232
Question: How should we select the correct kmer size for kmergenie analysis? I have read in wiki page that 'largest k-mer size to consider'. In that case I am looking for valuable comments from you people.

Answer: Smaller than the length of your reads, larger than 20. Attempt many different values and merge the assemblies in the end.
Question: Does anyone know where could I find a comprehensive list of human CNV regions? I am working with 1000G variants and want to filter out variants within CNV regions, to see whether some pattern change or not.

Thanks,
Federico

Answer: The 1000 Genomes Consortium recently published an integrated SV/CNV map for the 1000 Genomes phase 3 data:

http://www.nature.com/nature/journal/v526/n7571/full/nature15394.html

As a VCF file:

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/integrated_sv_map/ALL.wgs.integrated_sv_map_v2.20130502.svs.genotypes.vcf.gz

Alternatively, Decipher also has a population CNV map in tab-delimited format:

https://decipher.sanger.ac.uk/files/downloads/population_cnv.txt.gz
Question: I have a design of 6 conditions with three replicates each, with expression levels for ~51000 transcription products. I'm using Limma's voom transformation to make my data approximately Gaussian, which by ocular pat-down seems to be the case:

The mean-variance estimation is visible here - kinda overdispersed in the mid range but it doesn't look critically bad to me:

I fit the linear model thusly:

```
> labels
[1] A A A C C C G G G L L L T T T U U U
Levels: A C G L T U
design <- model.matrix(~labels)
fit2 <- eBayes( lmFit(voomCounts, design) )
```

However, this yields the following p-value distribution from the pairwise t-tests (looks like a misspecified model)

The `fit2$F.p.vals` has only 10-15 genes *above* 0.01, the rest are absolutely tiny; and when I plot some of the genes called for differential expression they look very unremarkable. What am I doing wrong?

Thanks so much :)

Answer: I do not know why are you doing voom counts

but if you want differentially expressed genes from 6 samples, lets say 3 replicates each of cases and controls

`df <-` "is your dataframe containing 6 coloumns (3 for each ) of normalized expression data"

then just go for

```
library(limma)

groups<-as.factor(c(rep("Cases",3),rep("Control",3)))
design<-model.matrix(~0+groups)
colnames(design)=levels(groups)
fit<-lmFit(df, design)
cont.matrix<-makeContrasts(Cases-Control,
levels=design)
fit2<-contrasts.fit(fit, cont.matrix)
ebfit<-eBayes(fit2)
topTable(ebfit, coef=1)
topTable(ebfit, number=Inf, p=0.05, adjust.method="none", coef=1)

## topTable would be your list of DEG
```
Question: Hiho

I want to create a rarefaction curve using R i.e. vegan {rearecurve/speccaccum}.

This worked out quite well (please see here):

![][1])

Anyways, as you can see in the figures the curves don't reach saturation. So I know that i have to kind of fit or simulate the missing data but this is exactly the problem. How can I produce a fit or something similar that shows me how large my sample has to be to reach saturation in discovered species?

Thanks for your time and help!

ps. code:

(raremax <- min(rowSums(t(species))))
Srare <- rarefy(t(species), raremax)
plot(specnumber(t(species)), Srare, xlab = "Observed No. of Species", ylab = "Rarefied No. of Species")
abline(0, 1)
rarecurve(t(species), step = 20, sample = raremax, col = "blue", cex = 0.6)

[1]: https://i.imgur.com/7gPo7gv.png

Answer: My quite naive suggestion is to estimate the richness with some estimator like Chao1 (should be in Vegan package), then extrapolate your richness curve to get 90% of estimated diversity and check for the sample size. See [iNEXT R package][1] for details on extrapolating rarefaction curves.

[1]: http://chao.stat.nthu.edu.tw/blog/software-download/inext-r-package/