Newbie question here, sorry. I am trying to generate a counts matrix for paired-end non-stranded data. I feel pretty good about the steps I have taken so far, but this is the first time I have used this particular type of data. The counts matrix I have generated has allele specific reads which I wasn't anticipating, and the counts seem vary off, especially between the two alleles. When generating a counts matrix with allele specific information, is there any special steps I need to consider? Also is there a way to get non-allele specific read counts?
I think a potential issue maybe with samtools or ht-seq. Feedback appreciated. I am willing to switch to bioconductor/R tools also. I don't have too much experience outside with RNA-seq outside of these programs though.
samtools sort -n accepted_hits.bam sample_sorted
samtools view sample_sorted.bam >sample.sam
htseq-count -s no --idattr=Parent sample.sam organism.gff > sample.counts
What would be helpful for me to know: Does it look like I'm on the right track? If yes, why are my allele reads so different? If no, what steps did I overlook, and do you have suggestions?
•
written
by
kpr •
1 answer
WARNING! Phased genome :)
As you have correctly noticed, C. albicans is one of few species which has a phased genome, i.e. sequence from both parental chromosomes (but keep in mind that in particular case of C. albicans we don't know if all A chromosomes are from one parental, and B chromosomes from another, here is the paper that generated this assembly is you are curious).
To my knowledge, today most of the software and pipelines (especially mappers and variant callers) are generally not well adapted in dealing with phased data. So if you are not really interested in allele-specific expression, its better to use reference genome and annotations from one parental (just subset fasta file and gff). If you are interested in allele-specific expression and will use phased genome, then you should expect to see a huge number of multimapped reads, because parental are very similar in DNA sequence. To avoid this at some extent you can restrict the number of allowed mismatches in settings of your mapper for example to 1 or 0 (in latter case you will miss the reads with sequencing errors though). In STAR mapper, these options are --outFilterMismatchNmax, --outFilterMismatchNoverLmax. More details are in manual.
And to speed up things, I would recommend you to use STAR mapper with option --quantMode GeneCounts, it will automatically generate count table of reads identical to HTSeq.
For reference genome taken form CGD, you can use this simple script:
import re
import sys
import os
from Bio import SeqIO
A_hapl = open("A_hapl.fasta", "w")
B_hapl = open("B_hapl.fasta", "w")
for seq_record in SeqIO.parse(sys.argv[1], "fasta"):
name = strseq_record.id)
seq = str(seq_record.seq)
if re.search(r'.*chr.A.*',name):
A_hapl.write("%s%s\n"%(">",name))
A_hapl.write("%s\n"%(seq))
elif re.search(r'.*chr.B.*', name):
B_hapl.write("%s%s\n"%(">",name))
B_hapl.write("%s\n"%(seq))
else:
A_hapl.write("%s%s\n"%(">",name))
A_hapl.write("%s\n"%(seq))
B_hapl.write("%s%s\n"%(">",name))
B_hapl.write("%s\n"%(seq))
A_hapl.close()
B_hapl.close()
It will generate files with A haplotype and B haplotype (A_hapl.fasta, B_hapl.fasta, respectively).
For gff file from CGD, you just can do like:
grep -e "A_C_albicans" -e "chrM_C_alb" C_albicans_SC5314_A22_current_features.gff > A_hapl.gff
Hope this helps
So basically you need to remap your reads to one of the hyplotypes.
IN R:
install.packages("seqinr")
library(seqinr)
# Read in fasta file
fileName = 'CA.fasta'
the.file = read.fasta(fileName)
M = grep("M_", names(the.file))
B_alleles = grep("B_", names(the.file))
A_alleles = grep("A_", names(the.file))
write.fasta(the.file, names(the.file)[sort(c(B_alleles, M))], "CA_A.fasta", open = "w")
write.fasta(the.file, names(the.file)[sort(c(A_alleles, M))], "CA_B.fasta", open = "w")
Log in to answer this question.
-
Phased Genome, Want reads From Both Alleles
written by kpr •
-
How to get counts data for allele specific expression analysis (ASE) from aligned bam files (RNA-se…
written by Naresh D J •
-
DESeq2 input for DE analysis
written by 1769mkc
-
Tophat2 Accepted hits Visualization and HT-seq Count not agreeing
written by Krish Srinivasan •
-
zeroes in ht-seq output
written by komal.yasmin •
-
Error in STAR's bams when running ht-seq
written by greener •
-
HTSeq-count wrong reads count
written by ariesliu324 •
-
Using HTSeq for paired-end files
written by Whoknows
-
Ht-Seq Read Count And Strand-Specificity
written by thecuriousbiologist
-
A Question About Samtools And Calculating Depth
written by Sameet •
Different how? Could you post an example?
What is the organism and are you mapping to a reference genome containing allele-specific information?
I have a wild, mutant1, and mutant2, 4 replicates each, below is a subset of the counts matrix that was generated. The rows are genes, and I think the A and B allele are marked at the end of the gene name. This is only a small subset, but there are two concerns. The first is that I would expect the alleles to be fairly similar, in some cases the reads between alleles differ by several hundred counts. The second concern is I will occasional see an extreme observation (such as C2_09060C) in which one sample's count differs by several thousand for one replicate but not the others. For example, for C2_09060C Wild, Mutant1, and Mutant2 had about the same numbers for all observations and replicates expect for one sample (Mutant2 below). There are quite a few occurrences of this.
Candida albicans, and I believe the reference genome contains allele specific information.
Indeed, it seems the reference genome assembly is diploid. I don't think HTSeq is the right tool for the job, I would search for allele-specif expression pipelines, or read some papers dealing with RNAseq analysis of Candida albicans.
If I use a haploid reference genome, will it fix the issue?
See my answer to your earlier post C Albicans Genome Count Matrix Total Genes: by selecting just one chromosome set, you may be leaving out a small number of features (and possibly of reference DNA, if there are structural variants between chromosome sets A and B). How small are the differences, I don't know.