cannin/biostars_qa huggingface dataset:
Question: Can I do differential gene expression analysis by using the cases and control data from two different GSEs?

Answer: If what you are asking is whether you can compare control and treatment across different experiments then the answer is that you have to more careful (while you are doing it) and more transparent (when disseminating the results) about what you are doing.

The statistical power of most methods will be greatly reduced and you would need to perform more careful error analysis and estimation.

What you are asking about is a variant of meta-analysis, where one combines the results of multiple, prior scientific studies to derive new insights.
Question: Hi,

I have multi-sample vcf and in this vcf, there are many sites which have uncalled or missing genotype. Is there a way to check which sample has greater number of uncalled genotypes in vcf. So that I can exclude that sample from further analysis.

Answer: A one liner using bioalcidaejdk: http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html

$ java -jar dist/bioalcidaejdk.jar -e 'stream().flatMap(G->G.getGenotypes().stream()).filter(G->!G.isCalled()).map(G->G.getSampleName()).collect(Collectors.groupingBy(Function.identity(), Collectors.counting())).forEach((K,V)->println(K+"\t"+V));' src/test/resources/test_vcf01.vcf | sort -t $'\t' -k2,2n

S3	8
S4	9
S5	14
S6	18
S2	23
S1	73

* `stream(). `get a stream of variants
* `flatMap(G->G.getGenotypes().stream()).` map to a stream of genotypes
* `filter(G->!G.isCalled()).` keep the uncalled genotype
* `map(G->G.getSampleName()). ` map to the sample name
* `collect(Collectors.groupingBy(Function.identity(), Collectors.counting())) ` convert to associative array sample/count
* .`forEach((K,V)->println(K+"\t"+V)); ` print the results.
Question: I have some RNAseq data - n=2, treated vs control. I have carried out alignment & quantification via Hisat2 > Stringtie > PrepDE.py > DESeq2

I would like to filter out D.E. genes that show high variance between the two replicates. But I don't seem to have an objective way of doing this.

![](https://image.ibb.co/hyyo3J/REP_CONCORDANCE.png)

In the screenshot as shown, I'd like to retain data that looks like the top panel, while filtering out data that looks like the bottom panel. As you can see, genes with very poor concordance between reps also show up with good p-vals and adjusted p-vals, and I was wondering why that is so.

Thanks

Answer: Check out the DESeq2 manual for Cooks Distance [here][1]. Cooks Distance is not applied to datasets with less than three replicates, but you can look at Mike Love's code sample [here][2] for an idea of how to apply it to your dataset.

maxCooks <- apply(assays(dds)[["cooks"]], 1, max)

From the Manual:

> The results function automatically flags genes which contain a Cook’s
> distance above a cutoff for samples which have 3 or more replicates.
> The p values and adjusted p values for these genes are set to NA. At
> least 3 replicates are required for flagging, as it is difficult to
> judge which sample might be an outlier with only 2 replicates. This
> filtering can be turned off with results(dds, cooksCutoff=FALSE)

[1]: https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
[2]: https://support.bioconductor.org/p/68277/
Question: Hi everyone.

I want to simplify the names of the samples of a VCF file by removing the suffix that was added during sequenciation but i am not being successful. Right now my samples are named like this:

123_4:A7TT6AKXX:1:1000000001 or ABC100:A7TT6AKXX:1:1000000001

In the first case I want to keep only the number before the underscore and in the second everything before the first colon (:). How can I do it?

Here's how the header is organized in my vcf file:

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 102_3:C6RL4ANXX:1:250529518 206E:C6RL4ANXX:2:250529583 AN07:C6RL4ANXX:2:250529585 B103:C6RL4ANXX:2:250529647 CA:C6RL4A

Thanks in advance.

Answer: bcftools **reheader** https://samtools.github.io/bcftools/bcftools.html

> -s, --samples FILE
>
> new sample names, one name per line, in the same order as they appear in the VCF file. Alternatively, only samples which need to be renamed can be listed as "old_name new_name\n" pairs separated by whitespaces, each on a separate line. If a sample name contains spaces, the spaces can be escaped using the backslash character, for example "Not\ a\ good\ sample\ name".
Question: I would like to measure the similarity among set of operons from different species. Is there a well-known method for this? If you can suggest me some good articles addressing this problem, I will be thankful to you.

Answer: Hi, please try [GeneGrouper][1]

It will look for all gene clusters in a set of genomes that contain a gene of interest. It will then show how similar those gene clusters are between taxa.

[1]: https://github.com/agmcfarland/GeneGrouper
Question: Hello,

I am trying to match Gencode's annotations to assemblies.

It is my understanding that the sequence of reference chromosomes changes only when there is a major version update (e.g. GRCh37 -> GRCh38). For minor versions (such as GRCh38.p2), patches (deltas between the major version and the new minor version) may be added (as well as haplotypes etc).

Gencode releases the following annotation: ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_22/gencode.v22.chr_patch_hapl_scaff.annotation.gtf.gz

that matches the following assembly: ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_22/GRCh38.p2.genome.fa.gz

If one doesn't want the patches, he can refer to the primary assembly: ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_22/GRCh38.primary_assembly.genome.fa.gz

which matches the following annotation: ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_22/gencode.v22.primary_assembly.annotation.gtf.gz

But then, what is this annotation for? ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_22/gencode.v22.annotation.gtf.gz

According to the description, this annotation describes reference chromosomes only. So why isn't this suitable for the primary assembly?

Also, it is often suggested not to mix and match Ensembl or Gencode annotations with UCSC assemblies, but given that there are 1:1 matchings (such as hg38 = GRCh38) that should be doable, as long as one takes care that chromosome names follow the same convention.

Similarly, if we remove patches and alternate loci from GRCh38.p2, at that point wouldn't we get back to the primary assembly GRCh38, except for differences in scaffolds? Then, if our annotations of choice only describe reference chromosomes, then those annotations, originally meant for GRCh38, would also work fine with GRCh38.p2. Isn't that the case?

Thank you for your help!

Roberto

Answer: The only thing that are added between patch versions (i.e., `.p1`, `.p2`) are patches. All of the haplotypes were present in GRCh38.0 and none will be added. New haplotype-like features are called "novel patches", until they are incorporated into a new genome assembly, GRCh39 - although that is likely to be a graph genome so the whole concept of haplotypes will be out the window. So removing patches from GRCh38.p2 would take you back to GRCh38.0, but the haplotypes have to stay. Therefore, the primary assembly for both will be identical, but the primary assembly of GRCh38.p2 is not the same as the complete genome for GRCh38.
Question: I wrote the below script which takes a SAM file as input and creates for each contig a SAM and FASTA file:

#!/usr/bin/env python3

import click
import os

@click.command()
@click.option('--sam', help="SAM files", required=True)
@click.option('--output', help="Output directory", required=True)
def retrieve_contig_reads(sam, output):

with open(sam) as align:
for line in align:
try:
parts = line.rstrip().split('\t')
illumina_id = parts[0]
Illumina_seq= parts[9]
contig_name = parts[2]

with open(os.path.join(output, contig_name + ".fasta"), 'a') as fasta:
fasta.write(">" + illumina_id + '\n')
fasta.write(Illumina_seq + '\n')

with open(os.path.join(output, contig_name + ".sam"), 'a') as sam_line:
sam_line.write(line + '\n')

except IndexError:
continue

if __name__ == '__main__':

retrieve_contig_reads()

Unfortunately, it ran for more than 2 weeks and has still not completed. Is there a faster way to retrieve the alignment and the reads for each contig?

Thank you in advance,

Answer: I wouldn't use python for this...

samtools index alignment.bam
samtools idxstats alignment.bam | cut -f 1 | grep -v '*' | parallel -j 8 --bar 'samtools view -o {}.bam alignment.bam {} && samtools fasta {}.bam | gzip > {}.fasta.gz'

1. Create an index of a your bam file
2. using idxstats, cut and grep -v get all contigs
3. use gnu parallel to parallelize processes, adjust `-j` accordingly to your system
4. for each contig: use samtools view to get the contig in a bam file, and samtools fasta to get the reads in fasta format
Question: Dear all,

I have paired-end fastq data generated with Illumina bcl2fastqv2.19 & sequenced on a Novaseq.The i5index is 7bp long, the i7 8bp long

R1.fastq.gz contains R1 101bp reads:

@A00154:125:HGKTMDMXX:1:1101:10420:1000 1:N:0:AACTGAGG+ATGCGTC
CTGGCCGTCTCAGCCGAGAAGCCGAGGATTGAATGGGCATGGAGACTGAACTACCCCTCTCACCTTTAGAGGTGGCTCCTCCAAGTCGGGGTTGACGCCCG
+
FFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

R2.fastq.gz contains 6bp UMI sequence

@A00154:125:HGKTMDMXX:1:1101:10420:1000 2:N:0:AACTGAGG+ATGCGTC
GCGCGT
+
FFFFFF

R3.fastq.gz contains R2 101bp reads:

@A00154:125:HGKTMDMXX:1:1101:10420:1000 3:N:0:AACTGAGG+ATGCGTC
CTTCATAGGCCACAAAAAGCCCATATATCAGTGTCATCCACTAAGCCTCAGACACTGCAGCACGGGCAGCGGCAGTGCCAGCTTCGCCCACACTGCCCCTC
+
FFFFFFFFFFFFFFFFFFFFFF:FF:FFF:FFFFFF:FFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

In a downstream analysis I want to use UMI-tools for deduplication. However for that I need the UMI be part of the read name.
@Instrument:RunID:FlowCellID:Lane:Tile:X:Y:**UMI** ReadNum:FilterFlag:0:IndexSequence or SampleNumber

There are tools to add a UMI to the read name when the UMI is present in the read itself. But in my case, the UMI is in a seperate fastq. How could this be achieved?

Answer: If you align the R1 and R3 to the genome of your choice, you can annotate it with the UMI using the fgbio tool AnnotateBamWithUmis: https://fulcrumgenomics.github.io/fgbio/tools/latest/AnnotateBamWithUmis.html, this does require you to have sufficient memory to store the entire R2 (UMI) .fastq file.
Question: Hi everyone,

This post is the second part of a question I already ask and I couldn't find any answers :

https://www.biostars.org/p/206213/#206286

I need a alternative to Cufflinks in order to process the bam files produced by Tophat2.
Indeed, I tried several time to use it, but each time, Cufflinks froze at the exact same location on my canu assembly (tig00000382:15043-15157) , and it can last forever.

My RNAseq data seem to be fine, since a colleague could analyse them (with an other assembly reference to align against with tophat), and use Cufflinks who didn't behave the same way, it didn't froze and do the job.

You can have more information about my problem on an issue I posted on github : https://github.com/cole-trapnell-lab/cufflinks/issues/72

It seems that I'm not the only one you got this problem.

For an alternative, I heard about Trinity, but I don't know that much about it, and about the differences between Trinity and Cufflinks.

Can I have some advices please ?

Thanks for helping !

Cheers,

Roxane

Answer: [StringTie][1] is pretty much a complete replacement for Cufflinks for transcript assembly and quantification. Its way faster than Cufflinks (e.g. a recent job that took Cufflinks 40 CPU hours took StringTie 1 CPU hour).

If you want to create your own transcript annotations, I definitely recommend using it.

For quantitaition, I'd recommend the alignment free algorithms (Kallisto, Salmon etc), as several people's benchmarks have recently shown them to be more accurate than the alignment based methods. See [here][2], [here][3] and [here][4].

[1]: https://ccb.jhu.edu/software/stringtie/
[2]: https://cgatoxford.wordpress.com/2016/08/17/why-you-should-stop-using-featurecounts-htseq-or-cufflinks2-and-start-using-kallisto-salmon-or-sailfish
[3]: http://f1000research.com/slides/4-741
[4]: http://www.ncbi.nlm.nih.gov/pubmed/26335491
Question: I read python basic codes and I am fimiliar with for loops, But not able to understand how to make a python programme to run hisat2 on multiple pair-end reads.

I want to run hisat2 mapping tool on mutiple pair end reads to map these on Genome.

I already have prepared genome index.

1)my input will be hisat2 location and genome index, pair-end read file
2) How should I give multiple hisat2 arguments and how to store sam file in a separate folder for each pair-end reads.

Can you help to make this python script ?

Location of hisat2:
/home/yog/software/hisat2-2.0.4/hisat2

Location of genome index:
data/Analysis/Radish_index

Location of pair end files:

/data/Analysis/mapped/

multiple pair end files like this:

216_5W_Ca1_R1.fq

216_5W_Ca1_R2.fq

218_5W_Ca1_R1.fq

218_5W_Ca1_R1.fq

216_5W_Co1_R1.fq

216_5W_Co1_R2.fq

218_5W_Co1_R1.fq

218_5W_Co1_R1.fq

Now I am running hisat2 on individual files by using this commonad:

/home/yog/software/hisat2-2.0.4/hisat2 -p 8 -x /data/Analysis/Radish_index --fr -1 /data/Analysis/mapped/216_5W_Ca1_R1.fq -2 /data/Analysis/mapped/216_5W_Ca1_R2.fq -S 216_5W_Ca1.sam

Answer: It would be easier to do this in bash:

for r1 in /data/Analysis/mapped/*R1.fq
do
#replace R1 by R2
r2=${r1/R1/R2}
# snip off the fq and add sam to the output file name
out=${r1%%_R1.fq}.sam
/home/yog/software/hisat2-2.0.4/hisat2 -p 8 -x /data/Analysis/Radish_index --fr -1 $r1 -2 $r2 -S $out
done

In Python it's similar:

import glob
import os
for r1 in glob.glob('/data/Analysis/mapped/*R1.fq'):
r2 = r1.replace('R1', 'R2')
out = r1.replace('_R1.fq', '.sam')
os.popen('/home/yog/software/hisat2-2.0.4/hisat2 -p 8 -x /data/Analysis/Radish_index --fr -1 %s -2 %s -S %s'%(r1, r2, out))

None of those are multithreaded - you can use the subprocess library instead of os.popen in Python so that each job runs in the background instead of it running each job one after the other, and you can end the command using & in bash to run all of the jobs at the same time, like this:

/home/yog/software/hisat2-2.0.4/hisat2 -p 8 -x /data/Analysis/Radish_index --fr -1 $r1 -2 $r2 -S $out &

Here I'm also not doing anything special to capture the standard error or the standard output, it would be better for you to have a look at all of the standard error.
Question: Hi,

I'm working with Omastandalone v. 2.3.0

I'm doing the all vs all comparisons with 25 species, they are almost done, however I got the next error recurrently, first just with one specie (then I deleted all the comparisons with that specie and they compute again right), but now they keep stopping in the last step with another species.

They are taking more time that was expected, and I can't figure it out how to solve this error.

Error, 'Matches between Echinops_telfairi/Balaena_mysticetus have been computed\
for different version of Echinops_telfairi'
/home/nis/daniela/OMA/bin/oma: line 196: kill: (24018) - No such process

Summary of OMA standalone All-vs-All computations:
--------------------------------------------------
Nr chunks started: 25 (0.02%)
Nr chunks finished: 132348 (99.92%)
Nr chunks finished w/o exported genomes: 132348 (99.92%)

I used Omastandalone before with another set of data and this error never came before.

I hope you could help me with this problem.

Thanks,

Answer: Hi,

the error is raised if the file in Cache/DB/Echinops_telfairi.db has been changed since it has been created (right in the beginning from the the fasta file). differences are based on computing a hash value from the file, so any modification of the file (for example changing an protein ID, etc) will cause this error to happen.

It's a bit strange that this happens so late when you are almost done with the All-vs-all computations. To me this really indicates that you somehow changed the file. If it is indeed a change on the file but not on any sequence, you can remove in the `Cache/AllAll/Echinops_telfairi/*` files the calles to `AssertDatabaseVersionsInSync`.

Best wishes
Adrian
Question: Hi,
I have an RNAseq data from TCR alpha and bed chain single cell sequencing which I am trying to analyze. I looked at the post https://www.biostars.org/p/148371/ and also read the Mark Davis paper and the protocol, which I thought would be appropriate to follow. But I cannot either find the dataset from the Mark Davis paper not the tool they used (vdjfasta). I was writing to ask what is the latest tool that would be appropriate for assembling the reads and getting a summary the TCR clones from my dataset.
Thanks,
- Pankaj

Answer: One of our post-docs has had nice results from [mixcr](https://github.com/milaboratory/mixcr). This was with an amplicon sequencing dataset, so I can't say how similar it is to your setup.
Question: Hi everyone,

I'm trying to analyze my counts data with DESeq2 and based on the tutorial of GSEA, DESeq2 has an output format that can be used directly in the GSEA ([here][1]). However, I'm reading their [workflow][2] and I don't find how to make this gct output format.
Any help is appreciated.

[1]: http://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/Using_RNA-seq_Datasets_with_GSEA
[2]: http://master.bioconductor.org/packages/release/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html#the-deseqdataset-object-sample-information-and-the-design-formula

Answer: Hi!

In this case you have to prepare the `gct` file "by hand" for analyze your data with GSEA. Once you have created your DESeq object (`dds`) you are going to retrieve the normalized counts based on the DESeq method:

norm_counts <- counts(dds, normalized = T)

Then, you must arrange this normalized counts matrix according to the minimum requirements for a [GCT](https://www.genepattern.org/file-formats-guide#GCT) file.
I suggest you to codificate the `norm_counts` as a `data.frame` and add an extra column called description (this column could contain the id's of your genes):

norm_counts <- as.data.frame(norm_counts)
norm_counts$description <- norm_counts$your_id_column

Move the "description" column to be the second in the df and save it using the `write.table` command:

write.table(norm_counts, "path_to_save_file.gct", sep = "\t", quote = F, row.names = F)

Finally, using a plain text editor put all the additional information (following the directions of the link) to have a gct file for GSEA.
For a detailed tutorial about `gct` and `cls` files watch [this video](https://www.youtube.com/watch?v=KY6SS4vRchY&t=2563s) but please follow the directions of the [link](https://www.genepattern.org/file-formats-guide#GCT)

Hope it could help!

Best regards!

Rodo
Question: I have 20 RNA-Seq samples obtained from normal cells (without any treatment). I want to identify from the transcript quantification data whether for a certain gene "A" which of the isoforms (A1, A2, A3, A4, A5) show differential expression or whether there is a differential transcript usage. The idea behind this comparison is to identify whether a regulatory element present in some isoforms of this gene (and not in other isoforms of the same gene) can affect the expression or usage of these transcripts.
In other words, I'm interested in finding out whether there is a significant difference in the expression of specific transcripts of the same gene in one condition.
I run "salmon" on the 20 fastq files to obtain the TPM of the transcripts (in my case there are ~1400 transcripts in ~500 genes).
However, I don't know how to perform this kind of comparisons or if there any R package that I could use to manage and compare TPM among transcripts of the same gene.
Any help will be appreciated and thanks in advance

Answer: This F1000 Article shows how to perform DTU analysis using Salmon quantification - [Swimming downstream: statistical analysis of differential transcript usage following Salmon quantification][1]

This makes use of classic packages such as DEXSeq, but optimised to use alignment-free quantifications.

Your hypothesis is pretty unclear, so I'm not sure what you're testing the different usage between. This article has the framework to make use of your input data, and perform the analysis you're interested in.

[1]: https://f1000research.com/articles/7-952
Question: Hi everyone.

I'm calling some variants using samtools (and bcftools) from a bwa aligned and sorted BAM.

For some reason, the samtools mpileup is reporting all zero quality scores, but I know the base and read quality scores in the BAM are good (viewed in IGV)

samtools mpileup -uvB -t DP -f ref.fa -r chrX:48,902,600-48,902,700 mapped_sorted.bam

gives me

```
CHROM POS ID REF ALT QUAL
chrX 48902600 . C 0
..
```

One for every base in the specified range

Here's a variant position:

chrX 48902688 . A C, 0

Base phred quality scores as reported in IGV range from 20-40 and read mapping quality is all 60.

Anyone run into this before?

Thanks!

Answer: I also spent last week or so trying to work through the project on calling variants and identifying causative gene (plants have been selected through phenotype. haven't finished the projects yet). I don't understand exactly what happens behind the scene, but `samtools mpileup` doesn't call variants! You need to use `bcftools call` to call variants. Here is some what helpful link http://samtools.sourceforge.net/mpileup.shtml Except be wary because it is out of date, bcftools have changes since that ti,me and you need to use `bcftools call` instead of `view`. Here is what I did to get my VCF files with quality scores

for i in Bams/*.bam; do (samtools mpileup -g -f reference.fa $i | bcftools call -c -v - > calledVCF/$(basename $i .bam).vcf)& done
Question: I have a large bam file that I need to subset/split into separate bam files by unique entries of one of the optional TAG:TYPE:VALUE fields. As a simplified example, say these are the optional fields in my example.bam file (i.e. fields following the mandatory 1-11 fields):

NH:i:6 HI:i:3 AS:i:94 nM:i:1 NM:i:1 CR:Z:TGCCCATCACTCAGGC CY:Z:CCCBBFGGGGGGGGGG CB:Z:TGCCCATCACTCAGGC-1 UR:Z:ACGCCCAGGT UY:Z:GGGGGGGGGG UB:Z:ACGCCCAGGT BC:Z:TGCGCAGC QT:Z:CCCCCGGG RG:Z:XXXXXX

NH:i:6 HI:i:3 AS:i:94 nM:i:1 NM:i:1 CR:Z:TGCCCATCACTCAGGC CY:Z:CCCBBFGGGGGGGGGG CB:Z:TGCCCATCACTCAGGC-1 UR:Z:ACGCCCAGGT

NH:i:7 HI:i:3 AS:i:96 nM:i:0 NM:i:0 CR:Z:TGCCCATCACTCAGGC CY:Z:BCBBCFGGGGGGGGGG CB:Z:TGCCCATCACTCAGGC-1 UR:Z:ACGCCCAGGT

NH:i:7 HI:i:3 AS:i:96 nM:i:0 NM:i:0 CR:Z:TGCCCATCACTCAGGC CY:Z:CCCCCGGGGGGGGGGG CB:Z:TGCCCATCACTCAGGC-1 UR:Z:ACGCCCAGGT

NH:i:3 HI:i:2 AS:i:96 nM:i:0 NM:i:0 CR:Z:TGCCCATCACTCAGGC CY:Z:CCCCCGGGGGGGGGGG CB:Z:TGCCCATCACTCAGGC-1 UR:Z:ACGCCCAGGT

NH:i:3 HI:i:2 AS:i:96 nM:i:0 NM:i:0 CR:Z:TGCCCATCACTCAGGC CY:Z:CCCCCGGGGGGGGGGG CB:Z:TGCCCATCACTCAGGC-1 UR:Z:ACGCCCAGGT

NH:i:6 HI:i:4 AS:i:96 nM:i:0 NM:i:0 CR:Z:CCGGTAGGTCATACTG CY:Z:CCCCCGGGGGGGGGGG CB:Z:CCGGTAGGTCATACTG-1 UR:Z:GCCGCCTTCT

NH:i:6 HI:i:4 AS:i:96 nM:i:0 NM:i:0 CR:Z:CCGGTAGGTCATACTG CY:Z:CCCCCGGGGGGGGGGG CB:Z:CCGGTAGGTCATACTG-1 UR:Z:GCCGCCTTCT

NH:i:6 HI:i:4 AS:i:96 nM:i:0 NM:i:0 CR:Z:CCGGTAGGTCATACTG CY:Z:CCCCCGGGGGGGGGGG CB:Z:CCGGTAGGTCATACTG-1 UR:Z:GCCGCCTTCT

NH:i:6 HI:i:4 AS:i:96 nM:i:0 NM:i:0 CR:Z:CCGGTAGGTCATACTG CY:Z:BBBCBGGEGGGGGGGG CB:Z:CCGGTAGGTCATACTG-1 UR:Z:GCCGCCTTCT

NH:i:6 HI:i:4 AS:i:96 nM:i:0 NM:i:0 CR:Z:CCGGTAGGTCATACTG CY:Z:CCCCCGGGGGGGGGGG CB:Z:CCGGTAGGTCATACTG-1 UR:Z:GCCGCCTTCT

How would I split this bam file by unique CB:Z: entries, such that the results would be 2 bam files as follows:

TGCCCATCACTCAGGC-1.bam (contains 6 reads)
CCGGTAGGTCATACTG-1.bam (contains 5 reads)

Thanks!
D

Answer: run a loop to find the distinct tag and pipe this into a samtools+awk to filter. Something like (not tested)

samtools view in.bam | cut -f 12- | tr "\t" "\n" | grep "^CR:Z:" | cut -d ':' -f 3 | sort | uniq | while read S; do samtools view -h in.bam | awk -v tag="CR:Z:$S" '($0 ~ /^@/ || index($0,tag)>0)' > ${S}.sam ; done