gt: 		for (CpGStatistics cpg : cpGStatisticsList) {
context: package edu.cwru.cbc.BSPAT.core;
  
  import edu.cwru.cbc.BSPAT.DataType.*;
  import org.apache.commons.lang3.StringUtils;
  import org.apache.commons.lang3.tuple.ImmutablePair;
  import org.apache.commons.lang3.tuple.Pair;
  import org.apache.commons.math.MathException;
  import org.apache.commons.math.distribution.NormalDistributionImpl;
  
  import java.io.IOException;
  import java.util.*;
  
  /**
   * Bisulfite sequences analysis. Include obtaining methylation string, methyl&SNP string. And group sequences by pattern.
   *
   * @author Ke
   */
  public class BSSeqAnalysis {
  	public static final int DEFAULT_TARGET_LENGTH = 70;
  	public static final double ASM_PATTERN_THRESHOLD = 0.2;
  	public static final double ASM_MIN_METHYL_DIFFERENCE = 0.2;
  	public static final double MEMU_PATTERN_THRESHOLD = 0.1;
  
  	/**
  	 * Execute analysis.
  	 *
  	 * @return List of ReportSummary
  	 */
  	public List<ReportSummary> execute(String experimentName,
  	                                   Constant constant) throws IOException, InterruptedException, MathException {
  		List<ReportSummary> reportSummaries = new ArrayList<>();
  		String inputFolder = constant.mappingResultPath + experimentName + "/";
  		String outputFolder = constant.patternResultPath + experimentName + "/";
  		ImportBismarkResult importBismarkResult = new ImportBismarkResult(constant.originalRefPath, inputFolder);
  		Map<String, String> referenceSeqs = importBismarkResult.getReferenceSeqs();
  		List<Sequence> sequencesList = importBismarkResult.getSequencesList();
  		if (sequencesList.size() == 0) {
  			throw new RuntimeException("mapping result is empty, please double check input!");
  		}
  
  		// 1. read refCoorMap / targetCoorMap
  		Map<String, Coordinate> refCoorMap = IO.readCoordinates(constant.coorFilePath, constant.coorFileName);
  		Map<String, Coordinate> targetCoorMap = getStringCoordinateMap(constant, referenceSeqs, refCoorMap);
  
  		// 2. group seqs by region
  		Map<String, List<Sequence>> sequenceGroupMap = groupSeqsByKey(sequencesList, Sequence::getRegion);
  
  		// 3. perform analysis for each region
  		for (String region : sequenceGroupMap.keySet()) {
  			ReportSummary reportSummary = new ReportSummary(region);
  			List<Sequence> seqGroup = sequenceGroupMap.get(region);
  			String refSeq = referenceSeqs.get(region);
  			Coordinate refCoor = refCoorMap.get(region);
  			Coordinate targetCoor = targetCoorMap.get(region);
  			if (!refCoor.getStrand().equals(targetCoor.getStrand())) {
  				throw new RuntimeException("target strand is not same to reference!");
  			}
  
  			// calculate target offset to reference
  			int targetStart, targetEnd;
  			switch (refCoor.getStrand()) {
  				case "+":
  					targetStart = targetCoor.getStart() - refCoor.getStart(); // 0-based
  					targetEnd = targetCoor.getEnd() - refCoor.getStart(); // 0-based
  					break;
  				case "-":
  					targetStart = refCoor.getEnd() - targetCoor.getEnd(); // 0-based
  					targetEnd = refCoor.getEnd() - targetCoor.getStart(); // 0-based
  					break;
  				default:
  					throw new RuntimeException("invalid reference strand!");
  			}
  			if (targetStart < 0 || targetEnd < 0 || targetStart >= targetEnd || targetStart >= refSeq.length() || targetEnd >= refSeq
  					.length()) {
  				throw new RuntimeException("invalid target coordinates! out of reference!");
  			}
  			int tmpStart = 0, tmpEnd = 0;
  			// used to include one more bp on both ends to detect CpG sites.
  			if (targetStart != 0) {
  				tmpStart--;
  			}
  			if (targetEnd != refSeq.length() - 1) {
  				tmpEnd++;
  			}
  			String tmpTargetRefSeq = refSeq.substring(targetStart + tmpStart, targetEnd + tmpEnd + 1);
  			int startCpGPos = tmpTargetRefSeq.indexOf("CG") + targetStart + tmpStart;
  			int endCpGPos = tmpTargetRefSeq.lastIndexOf("CG") + targetStart + tmpStart;
  			String targetRefSeq = refSeq.substring(targetStart, targetEnd + 1);
  			boolean isStartCpGPartial = (startCpGPos == targetStart - 1);
  			boolean isEndCpGPartial = (endCpGPos == targetEnd);
  			int cpgBoundedStart = isStartCpGPartial ? startCpGPos + 1 : startCpGPos;
  			int cpgBoundedEnd = isEndCpGPartial ? endCpGPos : endCpGPos + 1;
  			String cpgRefSeq = refSeq.substring(cpgBoundedStart, cpgBoundedEnd + 1);
  
  			// processing sequences
  			for (Sequence sequence : seqGroup) {
  				sequence.processSequence(refSeq);
  			}
  
  			// seqs in seqGroup got changed in updateTargetSequences().
  			Pair<List<Sequence>, List<Sequence>> filteredCuttingResultPair = updateTargetSequences(seqGroup,
  					targetStart, targetEnd, cpgBoundedStart, cpgBoundedEnd);
  			List<Sequence> CpGBoundSequenceList = filteredCuttingResultPair.getLeft();
  			List<Sequence> otherSequenceList = filteredCuttingResultPair.getRight();
  			reportSummary.setSeqOthers(otherSequenceList.size());
  
  			// sequence quality filter for CpGBoundSequenceList
  			reportSummary.setSeqCpGBounded(CpGBoundSequenceList.size());
  			Pair<List<Sequence>, List<Sequence>> filteredCpGSequencePair = filterSequences(CpGBoundSequenceList,
  					constant);
  			CpGBoundSequenceList = filteredCpGSequencePair.getLeft();
  			reportSummary.setSeqCpGAfterFilter(CpGBoundSequenceList.size());
  
  			// sequence quality filter for seqGroup
  			reportSummary.setSeqTargetBounded(seqGroup.size());
  			Pair<List<Sequence>, List<Sequence>> filteredTargetSequencePair = filterSequences(seqGroup, constant);
  			seqGroup = filteredTargetSequencePair.getLeft();
  			reportSummary.setSeqTargetAfterFilter(seqGroup.size());
  
  			List<Sequence> allMethylSequences = new ArrayList<>();
  			allMethylSequences.addAll(seqGroup);
  			allMethylSequences.addAll(CpGBoundSequenceList);
  
  			// if no sequence exist after filtering, return empty reportSummary
  			if (allMethylSequences.size() == 0) {
  				reportSummaries.add(reportSummary);
  				continue;
  			}
  
  			// calculate mismatch stat based on all sequences in reference region.
  			int[][] mismatchStat = calculateMismatchStat(targetRefSeq, targetStart, targetEnd, allMethylSequences);
  
  			// declare SNP
  			PotentialSNP potentialSNP = declareSNP(constant, allMethylSequences.size(),
  					mismatchStat, targetStart);
  
  			// write report
  			Report report = new Report(region, outputFolder, targetRefSeq, targetStart, cpgRefSeq, constant,
  					reportSummary);
  			report.writeReport(filteredTargetSequencePair, filteredCpGSequencePair, mismatchStat);
  
  			DrawPattern drawFigureLocal = new DrawPattern(constant.figureFormat, constant.refVersion,
  					constant.toolsPath, region, outputFolder, experimentName, targetCoorMap);
  
  			// generate methyl pattern output
  			List<Pattern> methylationPatternList = getMethylPattern(allMethylSequences, cpgBoundedStart,
  					cpgBoundedEnd);
  			System.out.println(experimentName + "\t" + region);
  			methylationPatternList = filterMethylationPatterns(methylationPatternList,
  					allMethylSequences.size(), StringUtils.countMatches(targetRefSeq, "CG"),
  					constant);
  			if (methylationPatternList.size() == 0) {
  				reportSummaries.add(reportSummary);
  				reportSummary.replacePath(Constant.DISKROOTPATH, constant.webRootPath, constant.host);
  				continue;
  			}
  			sortAndAssignPatternID(methylationPatternList);
  			reportSummary.addPatternLink(PatternLink.METHYLATION);
  			report.writePatterns(methylationPatternList, PatternLink.METHYLATION, allMethylSequences);
  			drawFigureLocal.drawPattern(reportSummary.getPatternLink(PatternLink.METHYLATION));
  
  			if (potentialSNP != null) {
  				// generate memu pattern output
  				List<Pattern> meMuPatternList = getMeMuPatern(seqGroup, methylationPatternList, potentialSNP,
  						targetStart, targetEnd);
  				meMuPatternList = filterPatternsByThreshold(meMuPatternList, seqGroup.size(),
  						MEMU_PATTERN_THRESHOLD);
  				if (meMuPatternList.size() != 0) {
  					reportSummary.addPatternLink(PatternLink.METHYLATIONWITHSNP);
  					report.writePatterns(meMuPatternList, PatternLink.METHYLATIONWITHSNP, seqGroup);
  					drawFigureLocal.drawPattern(reportSummary.getPatternLink(PatternLink.METHYLATIONWITHSNP));
  				}
  
  				// ASM
  				Pair<Pattern, Pattern> allelePatterns = getAllelePatterns(seqGroup, potentialSNP);
  				Pattern allelePattern = allelePatterns.getLeft();
  				Pattern nonAllelePattern = allelePatterns.getRight();
  				if (allelePattern.getSequenceMap().size() != 0 && nonAllelePattern.getSequenceMap().size() != 0) {
  					PatternResult patternWithAllele = patternToPatternResult(allelePattern, report.getCpgStatList(),
  							seqGroup.size(), targetStart);
  					PatternResult patternWithoutAllele = patternToPatternResult(nonAllelePattern,
  							report.getCpgStatList(), seqGroup.size(), targetStart);
  
  					if (!hasASM(patternWithAllele, patternWithoutAllele)) {
  						reportSummary.setHasASM(false);
  					} else {
  						reportSummary.setHasASM(true);
  						drawFigureLocal.drawASMPattern(reportSummary, patternWithAllele, patternWithoutAllele,
  								constant.logPath);
  					}
  				}
  			}
  
  			reportSummary.replacePath(Constant.DISKROOTPATH, constant.webRootPath, constant.host);
  			reportSummaries.add(reportSummary);
  		}
  		return reportSummaries;
  	}
  
  	private synchronized Map<String, Coordinate> getStringCoordinateMap(Constant constant, Map<String, String> referenceSeqs,
  	                                                                    Map<String, Coordinate> refCoorMap) throws
  			IOException {
  		Map<String, Coordinate> targetCoorMap = IO.readCoordinates(constant.targetPath, constant.targetFileName);
  		for (String key : refCoorMap.keySet()) {
  			// if no given targetCoor, get position of first CpG and generate DEFAULT_TARGET_LENGTH bp region.
  			if (!targetCoorMap.containsKey(key)) {
  				String refString = referenceSeqs.get(key);
  				String strand = refCoorMap.get(key).getStrand();
  				if (strand.equals("+")) {
  					int firstPos = refString.indexOf("CG");
  					if (firstPos == -1) {
  						throw new RuntimeException("no CG found in reference!");
  					} else {
  						// from first CpG to min(ref end, fisrt CpG + DEFAULT_TARGET_LENGTH)
  						int startPos = refCoorMap.get(key).getStart() + firstPos;
  						int endPos = refCoorMap.get(key).getStart() + firstPos + DEFAULT_TARGET_LENGTH - 1;
  						targetCoorMap.put(key,
  								new Coordinate(refCoorMap.get(key).getId(), refCoorMap.get(key).getChr(), strand,
  										startPos,
  										endPos < refCoorMap.get(key).getEnd() ? endPos : refCoorMap.get(key).getEnd()));
  					}
  				} else if (strand.equals("-")) {
  					int firstPos = refString.lastIndexOf("CG") + 2;
  					if (firstPos == -1) {
  						throw new RuntimeException("no CG found in reference!");
  					} else {
  						int startPos = refCoorMap.get(key).getEnd() - firstPos + 1;
  						int endPos = startPos + DEFAULT_TARGET_LENGTH - 1;
  						targetCoorMap.put(key,
  								new Coordinate(refCoorMap.get(key).getId(), refCoorMap.get(key).getChr(), strand,
  										startPos,
  										endPos < refCoorMap.get(key).getEnd() ? endPos : refCoorMap.get(key).getEnd()));
  					}
  				}
  			}
  		}
  		if (constant.getTargetFileName() == null) {
  			constant.targetFileName = "defaultTarget.coor";
  			IO.writeCoordinates(constant.targetPath + constant.targetFileName, targetCoorMap);
  		}
  		return targetCoorMap;
  	}
  
  	private PotentialSNP declareSNP(Constant constant, int totalTargetSeqenceCount, int[][] mismatchStat,
  	                                int targetStart) {
  		List<PotentialSNP> potentialSNPList = new ArrayList<>();
  		double threshold = totalTargetSeqenceCount * constant.SNPThreshold;
  		for (int i = 0; i < mismatchStat.length; i++) {
  			int count = 0;
  			if (mismatchStat[i][0] >= threshold) {
  				potentialSNPList.add(new PotentialSNP(i + targetStart, 'A'));
  				count++;
  			}
  			if (mismatchStat[i][1] >= threshold) {
  				potentialSNPList.add(new PotentialSNP(i + targetStart, 'C'));
  				count++;
  			}
  			if (mismatchStat[i][2] >= threshold) {
  				potentialSNPList.add(new PotentialSNP(i + targetStart, 'G'));
  				count++;
  			}
  			if (mismatchStat[i][3] >= threshold) {
  				potentialSNPList.add(new PotentialSNP(i + targetStart, 'T'));
  				count++;
  			}
  			if (mismatchStat[i][4] >= threshold) {
  				potentialSNPList.add(new PotentialSNP(i + targetStart, 'N'));
  				count++;
  			}
  			if (count >= 2) {
  				//more than two SNP allele in same position!
  				return null;
  			}
  		}
  		switch (potentialSNPList.size()) {
  			case 0:
  				return null;
  			case 1:
  				return potentialSNPList.get(0);
  			default:
  				// More than 1 SNP in the region!
  				return null;
  		}
  	}
  
  	private boolean hasASM(PatternResult patternWithAllele, PatternResult patternWithoutAllele) {
  		// use 0.2 as threshold to filter out unequal patterns. ASM pattern should be roughly equal.
  		if (patternWithAllele.getPercent() < ASM_PATTERN_THRESHOLD ||
  				patternWithoutAllele.getPercent() < ASM_PATTERN_THRESHOLD) {
  			return false;
  		}
  		List<CpGSitePattern> cgListWithAllele = patternWithAllele.getCpGList();
  		List<CpGSitePattern> cgListWithoutAllele = patternWithoutAllele.getCpGList();
  		// if there is at least one cpg site with different methyl type and the different bigger than 0.2, it is ASM
  		for (int i = 0; i < cgListWithAllele.size(); i++) {
  			if (cgListWithAllele.get(i).getMethylType() != cgListWithoutAllele.get(i).getMethylType() &&
  					Math.abs(cgListWithAllele.get(i).getMethylLevel() - cgListWithoutAllele.get(i).getMethylLevel()) >=
  							ASM_MIN_METHYL_DIFFERENCE) {
  				return true;
  			}
  		}
  		return false;
  	}
  
  	private PatternResult patternToPatternResult(Pattern pattern, List<CpGStatistics> cpGStatisticsList,
  	                                             int totalCount, int targetStart) {
  		PatternResult patternResult = new PatternResult();
  		Map<Integer, CpGSitePattern> cpGSiteMap = new HashMap<>();
  