public void createSAM(List<PeptideSearchResult> peptideSearchResults,
Writer output, Writer bedWriter)
throws PeptideSequenceGeneratorException, IOException, FastaParserException {
LOG.debug("creating sam file");
List<SAMEntry> samEntries = new ArrayList<SAMEntry>();
PeptideSequenceGenerator sequenceGenerator = new PeptideSequenceGeneratorImpl(
genome, proteinToOLNMap, chromosomeDir);
Set<String> foundProteins = new HashSet<String>();
PeptideSearchResultFilter peptideFilter = null;
if (confidenceScore != null) {
peptideFilter = new ConfidenceScoreFilter(confidenceScore);
}
for (PeptideSearchResult result : peptideSearchResults) {
if (peptideFilter != null && !peptideFilter.accepts(result)) {
continue;
}
PeptideSequence peptide = sequenceGenerator.getPeptideSequence(result);
if (peptide == null) {
LOG.warn("Error while geting peptide sequnce for " + result.getId());
continue;
}
String proteinName = result.getProteinName();
String resultName = proteinName + "." + result.getId();
int peptideStart = peptide.getStartIndex()
+ peptide.getGeneInfo().getStart();
if (bedWriter != null && !foundProteins.contains(proteinName)) {
foundProteins.add(proteinName);
BedLineOutputter bedLineOutputter = new BedLineOutputter(
peptide, proteinName);
bedWriter.write(bedLineOutputter.toString());
}
if (DebuggingFlag.get_use_mascot_score_flag() == 1) {
int mapq = result.getConfidenceScore()
.round(new MathContext(0)).intValue();
samEntries.add(new SAMEntry(resultName, peptide.getGeneInfo(),
peptideStart, peptide.getCigarString(), peptide
.getNucleotideSequence(), mapq));
} else {
int mapq = result.getConfidenceScore().round(new MathContext(0)).intValue();
String sequnece = peptide.getNucleotideSequence();
String outputSequence = (peptide.getGeneInfo().getDirection() == -1) ? new StringBuilder(StringUtils.replaceChars(sequnece, "ACGT", "TGCA")).reverse().toString() : sequnece;
SAMEntry entry = new SAMEntry(resultName, peptide.getGeneInfo(), peptideStart, peptide.getCigarString(), outputSequence, mapq);
entry.setChromosomeLength(sequenceGenerator.getFastaParser().getChromosomeLength(peptide.getGeneInfo().getChromosome()));
samEntries.add(entry);
}
try {
if (DebuggingFlag.get_sbi_debug_flag() == 1) {
CodonTranslationTable translationTable = CodonTranslationTable
.parseTableFile(translationTableFile);
String nucleotideString = peptide.getNucleotideSequence();
int direction = peptide.getGeneInfo().getDirection();
String mascotPeptideString = result.getPeptideSequence();
String predictedAminoAcidSequence = new String("");
if (direction != 1) {
StringBuilder invertedReversedSequence = new StringBuilder(
StringUtils.replaceChars(nucleotideString,
"ACGT", "TGCA")).reverse();
predictedAminoAcidSequence = translationTable
.proteinToAminoAcidSequence(invertedReversedSequence
.toString());
} else {
predictedAminoAcidSequence = translationTable
.proteinToAminoAcidSequence(nucleotideString);
}
if (!predictedAminoAcidSequence.equals(mascotPeptideString)) {
String samEntryStrng = new SAMEntry(resultName,
peptide.getGeneInfo(), peptideStart,
peptide.getCigarString(),
peptide.getNucleotideSequence()).toString();
LOG.info("Incorrect nucleotide sequence for following SAM entry:\n"
+ samEntryStrng);
}
}
} catch (Exception npe) {
LOG.info("Problem with internal validation of output nucleotide sequence!");
npe.printStackTrace();
System.err.println(npe);
}
}
String prevChromosome = null;
Collections.sort(samEntries, new SAMEntryComparator());
// ////////////////////////////////////////////////
// // Location to store SAM file headers ///
// ////////////////////////////////////////////////
//reverse if needed
//writing header
ArrayList<String> chromosomes = new ArrayList<String>();
for (SAMEntry samEntry : samEntries) {
if (chromosomes.size() > 0 && chromosomes.get(chromosomes.size() - 1).equalsIgnoreCase(samEntry.getRname())) {
continue;
}
chromosomes.add(samEntry.getRname());
}
output.write("@HD\tVN:1.0\n");
for (String chromosome : chromosomes) {
output.write("@SQ\tSN:" + chromosome + "\tLN:" + ((sequenceGenerator.getFastaParser().getChromosomeLength(chromosome))) + "\n");
}
for (SAMEntry samEntry : samEntries) {
String chromosome = samEntry.getRname();
if (!chromosome.equals(prevChromosome)) {