实验记录 | 关于sam文件的一些疑问解答

1。同一个序列名，为何比对的sequence会不一样？

(base) [xxzhang@mu02 trash]$ awk ' $2=="0" && (80029677< $4 )&& ($4<80029976 )' AluYe2_chr5.sam

A00928:207:HYLCHDSXY:2:1144:25898:22075 0       chr5    80029767        60      151M    *       0       0       GGCTAACATGGTGAAACCCCGTCTCTACTAAAAATACAAAAAAAAAAAAAAAAATTAGCCGGGCATGGTGGTGGGCACCTTTAGTCCCAGCTACTCGGGAGGCTGAGATAGGAGAATGGCGTGAACCCGGGAGGCGGAGGTTGCAGTAAGC      *       NM:i:0  MD:Z:151        AS:i:151        XS:i:97
A00928:207:HYLCHDSXY:2:1144:25898:22075 0       chr5    80029821        50      151M    *       0       0       TTAGCCGGGCATGGTGGTGGGCACCTTTAGTCCCAGCTACTCGGGAGGCTGAGATAGGAGAATGGCGTGAACCCGGGAGGCGGAGGTTGCAGTAAGCCGAGATCGCGCTGCTGCACTCCAGCCTGGGCGACAGCGAGACTCCATCTCAAAA      *       NM:i:0  MD:Z:151        AS:i:151        XS:i:128        XA:Z:chr8,-47957525,18M2D133M,5;
A00928:207:HYLCHDSXY:2:2136:5195:35603  0       chr5    80029771        6       151M    *       0       0       AACATGGTGAAACCCCGTCTCTACTAAAAATACACACACACAAAAAAAAATTAGCCGGGCATGGTGGTGGGCACCTTTAGTCCCAGCTACTCGGGAGGCTGAGATAGGAGAATGGCGTGAACCCGGGAGGCGGAGGTTGCAGTGAGCCGAG      *       NM:i:5  MD:Z:34A1A1A1A102A7     AS:i:126        XS:i:121        XA:Z:chr8,-47957577,110M7I34M,9;
A00928:207:HYLCHDSXY:2:2136:5195:35603  0       chr5    80029810        25      151M    *       0       0       ACAAAAAAAAATTAGCCGGGCATGGTGGTGGGCACCTTTAGTCCCAGCTACTCGGGAGGCTGAGATAGGAGAATGGCGTGAACCCGGGAGGCGGAGGTTGCAGTGAGCCGAGATCGCGCTGCTGCACTCCAGCCTGGGCGACAGCGAGACT      *       NM:i:2  MD:Z:1A102A46   AS:i:144        XS:i:132        XA:Z:chr8,-47957536,7M2D142M2S,4;

比如对于序列A00928:207:HYLCHDSXY:2:1144:25898:22075，它后面的序列完全不一样。为什么会出现这种情况，以及如何解释它？

(base) [xxzhang@mu02 trash]$ grep "A00928:207:HYLCHDSXY:2:1144:25898:22075" result.sam 
                                                                                                                                 
A00928:207:HYLCHDSXY:2:1144:25898:22075 99      chr5    80029767        60      151M    =       80029821        205     GGCTAACATGGTGAAACCCCGTCTCTACTAAAAATACAAAAAAAAAAAAAAAAATTAGCCGGGCATGGTGGTGGGCACCTTTAGTCCCAGCTACTCGGGAGGCTGAGATAGGAGAATGGCGTGAACCCGGGAGGCGGAGGTTGCAGTAAGC      FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:F::,FFFF:::FF::FFFFFFFFF,FFF,,::F,:F,::,,:FF::F::F:,,F::F:,:F:FFFFF:FF,FFFFF,:F::F::FFFFF::::FF,F NM:i:0  MD:Z:151        MC:Z:151M       AS:i:151     XS:i:97
A00928:207:HYLCHDSXY:2:1144:25898:22075 147     chr5    80029821        60      151M    =       80029767        -205    TTAGCCGGGCATGGTGGTGGGCACCTTTAGTCCCAGCTACTCGGGAGGCTGAGATAGGAGAATGGCGTGAACCCGGGAGGCGGAGGTTGCAGTAAGCCGAGATCGCGCTGCTGCACTCCAGCCTGGGCGACAGCGAGACTCCATCTCAAAA      FFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF NM:i:0  MD:Z:151        MC:Z:151M       AS:i:151     XS:i:93


(base) [xxzhang@mu02 family]$ grep "A00928:207:HYLCHDSXY:2:1144:25898:22075" redup.F2.sam

A00928:207:HYLCHDSXY:2:1144:25898:22075 99      AluYe2  90      11      37M11I103M      =       133     191     GGCTAACATGGTGAAACCCCGTCTCTACTAAAAATACAAAAAAAAAAAAAAAAATTAGCCGGGCATGGTGGTGGGCACCTTTAGTCCCAGCTACTCGGGAGGCTGAGATAGGAGAATGGCGTGAACCCGGGAGGCGGAGGTTGCAGTAAGC      FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:F::,FFFF:::FF::FFFFFFFFF,FFF,,::F,:F,::,,:FF::F::F:,,F::F:,:F:FFFFF:FF,FFFFF,:F::F::FFFFF::::FF,F AS:i:174        XS:i:181        XN:i:0  XM:i:10 XO:i:1       XG:i:11 NM:i:21 MD:Z:8C44G0A5C4G3G26G0C30C7G3   YS:i:212        YT:Z:CP MQ:i:11 MC:Z:148M3S     ms:i:5551
A00928:207:HYLCHDSXY:2:1144:25898:22075 147     AluYe2  133     11      148M3S  =       90      -191    TTAGCCGGGCATGGTGGTGGGCACCTTTAGTCCCAGCTACTCGGGAGGCTGAGATAGGAGAATGGCGTGAACCCGGGAGGCGGAGGTTGCAGTAAGCCGAGATCGCGCTGCTGCACTCCAGCCTGGGCGACAGCGAGACTCCATCTCAAAA      FFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF AS:i:212        XS:i:201        XN:i:0  XM:i:12 XO:i:0  XG:i:0       NM:i:12 MD:Z:10G0A5C4G3G26G0C30C7G14C0A32G5     YS:i:174        YT:Z:CP MQ:i:11 MC:Z:37M11I103M ms:i:4661

因为之前乱七八糟的处理方式，现在让我弄得都有些混乱了。有点爆炸。现在需要的就是把心情平复下来，不吃饭了，今天。
对于这个问题，其实就是这两条具有相同的序列名的sequence，就是我们说的双端测序的两条reads。
下面的两条信息，其实是双端测序的结果。而上面的两条信息，就是单端比对的结果。我们来解释其为比较的结果。这样看来，这个问题提的太傻了。
我现在好痛苦啊，一点食欲都没有。就感觉到情绪一直在下沉，一点快乐的想法都没有。

2。mapping 模式都是151M，为什么mapping的质量却有某种程度的差异？M（包括错配）和S的区别是什么？

A00928:207:HYLCHDSXY:2:1101:3088:1047   99      chr6    169751668       0       151M    =       169751835       318     GNTGGGTGTAGGGCCCGGTTCGATCCCGAGCTAGGCAGGGAGTCGGCGCCAGGCTGGGTGGTGCTCGGCTACGCGGAGTGGGCGAGCGAGCACGCGCCTGCGGTGGCCGGCGGTCCCGCGCTGGAGGGCGCTGGCGACGTCGCGGCCCTGG      F#,:FF,FFFF:F:FFFF::FFFFFFFFFFFFFFFFF:FFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFF:FFFF,FFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFF NM:i:2  MD:Z:1G0C148    MC:Z:151M       AS:i:148     XS:i:148
A00928:207:HYLCHDSXY:2:1101:3088:1047   147     chr6    169751835       0       151M    =       169751668       -318    GGACGCTCGGCCCTGCAAGACGCCTGGCGGGCCCTGCCGGCCTCCCTGGGCCAGGCTCGGTTTCTCCGTCGCCTCCTGGCCCTGAGCTGGAACGCGGCGGACGGTCAGGGAGGGCTGTGCGCGGTGGCGCGGGGGCGGAAAGCCGCAGGTC      FFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF NM:i:0  MD:Z:151        MC:Z:151M       AS:i:151     XS:i:151

mapping quality=0代表的意思是什么？

MAPQ: MAPping Quality. It equals −10 log10 Pr{mapping position is wrong}, rounded to the nearest
integer. A value 255 indicates that the mapping quality is not available.
这个mapping quality 的意思是说没有匹配上，匹配的质量不高是吗？就是说匹配的错误率为1。所以是否应该根据这个指标来过度呢？
而且虽然显示的是151M，但是这里其实是包括了错配的情况的。仔细想想插入和缺失确实也会出现，但是错配的情况的确用M来表示了。错配的程度和mapping quality的值怎样的相互影响的程度，需要仔细考虑一下。


（通过上图可以理解S、H、M的区别）

3。关于次优的比对结果的呈现方式。证明：比对到多个位置的reads并没有被过滤掉，而是结合其双端的比对情况，选择的是最优的比对的结果进行后续的分析。

按照我目前的了解，如果一条reads能够匹配到多个地方，那么在sam文件中会出现的方式是：
（1）同一个序列名，会有多条的比对的结果。（对于双端的reads，数量大于等于3）
（2）对于次优的比对结果，flag会大于256。

#比对情况1
#序列1
A00928:207:HYLCHDSXY:2:1101:14552:1047  99      chr19   56392542        60      51M100S =       56392542        51      GNCCACAAGTAGCCAAGAGGTTGAATCATTTTTGGCACAGCCATATAATGGCTGTCTCTTATACACATCTCCGAGCCCACGAGACCCGAGGCAATCTCGGTTTCCGTCTTCGGCGTGGAAAAGGGGGGGGGGGGGGGGGGGGGGGGGGGGG      F#FFFFFFFFFFFFFFFFFF::FFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,FF:FF,,F,,F,,,,,,,F,,,,,F:,,,,:FFFFFFFFFFFFFFFFFFFFFFFFFF NM:i:1  MD:Z:1T49       MC:Z:100S51M    AS:i:49      XS:i:20 SA:Z:chr14,77269795,-,30M121S,0,0;
#序列2
A00928:207:HYLCHDSXY:2:1101:14552:1047  147     chr19   56392542        60      100S51M =       56392542        -51     CCCCCCCCCTTTTTAATGATACGGCGACCACCGAGATCTACACACATTGCGTTTCCAGGCGCGTCTGTCGTCGGCAGCGTCAGATGTGTATAAGAGACAGGTCCACAAGTAGCCAAGAGGTTGAATCATTTTTGGCACAGCCATATAATGG      FFFFFFFF,,:FFFFFFFF:,FF:FFFFFFFFFFFFFFFFFFF:F:,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF NM:i:0  MD:Z:51 MC:Z:51M100S    AS:i:51 XS:i:20
#比对情况2
A00928:207:HYLCHDSXY:2:1101:14552:1047  2163    chr14   77269795        0       30M121H chr19   56392542        0       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCT  FFFFFFFFFFFFFFFFFFFFFFFFFF:,,,  NM:i:0  MD:Z:30 MC:Z:100H51M    AS:i:30 XS:i:30 SA:Z:chr19,56392542,+,51M100S,60,1;
#比对情况3
A00928:207:HYLCHDSXY:2:1101:14552:1047  355     chr2    32916441        0       121H30M chr19   56392542        0       *       *       NM:i:0  MD:Z:30 MC:Z:100H51M    AS:i:30
A00928:207:HYLCHDSXY:2:1101:14552:1047  355     chr2    32916486        0       121H30M chr19   56392542        0       *       *       NM:i:0  MD:Z:30 MC:Z:100H51M    AS:i:30
A00928:207:HYLCHDSXY:2:1101:14552:1047  355     chr2    32916351        0       121H30M chr19   56392542        0       *       *       NM:i:0  MD:Z:30 MC:Z:100H51M    AS:i:30
A00928:207:HYLCHDSXY:2:1101:14552:1047  355     chr2    32916401        0       121H30M chr19   56392542        0       *       *       NM:i:0  MD:Z:30 MC:Z:100H51M    AS:i:30

#比对情况1是最好的一种比对的情况。序列1和序列2都比对到同一个染色体上（按照预期或者采样的技术应该是这种情况的。），且比对的质量值很高。
#比对情况2是剪掉了一部分，然后和我们的参考基因组比对上了。
#比对情况3就是次有的比对结果们！但是这些情况，他们的质量值都很低，且flag>256,将会在我们之后的处理过程中被筛选掉。

所以，这样的一个尝试，足以可以证明：比对到多个位置的reads对保留到我们的结果文件中，我们在后续的处理的过程中，会把一些次优的结果过滤掉。然后保留最优的双端匹配的结果用于后续的分析。

4。双端序列中，比对到负链的那一条链，是反向互补，反向还是没有序列上的调整变化？即负链的reads在结果文件中是怎样的呈现形式。

还是以我们刚刚举例的那一条reads为例：A00928:207:HYLCHDSXY:2:1101:14552:1047

grep "A00928:207:HYLCHDSXY:2:1101:14552:1047" ATAC_S1_20210413NA_S2_L002_R1_001.fastq
@A00928:207:HYLCHDSXY:2:1101:14552:1047 1:N:0:CCGAGGCA
GNCCACAAGTAGCCAAGAGGTTGAATCATTTTTGGCACAGCCATATAATGGCTGTCTCTTATACACATCTCCGAGCCCACGAGACCCGAGGCAATCTCGGTTTCCGTCTTCGGCGTGGAAAAGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
+
F#FFFFFFFFFFFFFFFFFF::FFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,FF:FF,,F,,F,,,,,,,F,,,,,F:,,,,:FFFFFFFFFFFFFFFFFFFFFFFFFF
###########################################
grep "A00928:207:HYLCHDSXY:2:1101:14552:1047" -A 4 ATAC_S1_20210413NA_S2_L002_R3_001.fastq
@A00928:207:HYLCHDSXY:2:1101:14552:1047 3:N:0:CCGAGGCA
CCATTATATGGCTGTGCCAAAAATGATTCAACCTCTTGGCTACTTGTGGACCTGTCTCTTATACACATCTGACGCTGCCGACGACAGACGCGCCTGGAAACGCAATGTGTGTAGATCTCGGTGGTCGCCGTATCATTAAAAAGGGGGGGGG
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,:F:FFFFFFFFFFFFFFFFFFF:FF,:FFFFFFFF:,,FFFFFFFF

有趣的事情出现了。这个是我之前没有注意到的！我发现比对到负链的reads，在我们的sam文件中比逆反的排列了。也就是说，现在在sam文件中呈现的序列是原来的序列的反向互补的序列。也就是那条负链在正链上所对应的reads的序列的情况。

5。关于如何解释，R1比对正链，R2比对到正链这两种情况?以及R1和R2这两条reads之间的关系？

从这张示意图上看，R1测的是上面那一条链。R2测的是下面那一条链。不管怎样说R1和R2是配对出现的。肯定一条是正链，一条是负链。那么可能存在这种情况嘛？就是在R1中出现的reads 既可以比对到正链，又可以比对到负链？R2同理。因为在我们的结果中，我的确遇到了这种情况。

A00928:207:HYLCHDSXY:2:2410:21007:35211	99	chr5	15289	51	75M76S	=	15289	74	ACCACACAGGAGGTAGAGGTTGCAGTGAGCTGAGATTGCACCTTTGCACTCCAGCCTGGGCGACAGAGCCAGACCTGTCTCTTATACACATCTCCGAGCCCACGAGACCCGAGGCAATCGCGTATTCCCTCGTATGCTTGAAAAAGGGGGG	FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFF,FF:FFFFFFFFFFF:FFFFFFFFFF,:F,FF:FF,,,,,,,:,,,,,,,::F:,,,,,FF	NM:i:0	MD:Z:75	MC:Z:77S74M	AS:i:75	XS:i:69	XA:Z:chr1,+349618,6S69M76S,0;chr2,+113615216,6S69M76S,0;chr5,-181470707,76S69M6S,0;
A00928:207:HYLCHDSXY:2:2512:29270:19836	163	chr5	15284	60	80M71S	=	15284	79	GGTTTACCACACAGGAGGTAGAGGTTGCAGTGAGCTGAGATTGCACCTTTGCACTCCAGCCTGGGCGACAGAGCCAGACCTGTCTCTTATACACATCTGACGCTGCCGACGACAGACGCGGGTAGCTGAAGGATTTGTGTAGATCTCGGTG	FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF	NM:i:0	MD:Z:80	MC:Z:72S79M	AS:i:80	XS:i:69	XA:Z:chr1,+349618,11S69M71S,0;chr2,+113615216,11S69M71S,0;chr5,-181470707,71S69M11S,0;
A00928:207:HYLCHDSXY:2:2525:17047:22874	163	chr5	15250	60	151M	=	15543	444	GATTACAGGGTTCACAGCCTTAAAGGTTAGTAGAGGTTTACCACACAGGAGGTAGAGGTTGCAGTGAGCTGAGATTGCACCTTTGCACTCCAGCCTGGGCGACAGAGCCAGACCCTGTCTCAAAAAAAAAATTTTTTAAAGGAAAACTATA	FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFF:FFFFFFFFFFF:,FFFFFF:FF:FFFFFFFFFFFFFFF	NM:i:1	MD:Z:130T20	MC:Z:151M	AS:i:146	XS:i:93
A00928:207:HYLCHDSXY:2:2560:15881:15092	99	chr5	15289	51	75M76S	=	15289	74	ACCACACAGGAGGTAGAGGTTGCAGTGAGCTGAGATTGCACCTTTGCACTCCAGCCTGGGCGACAGAGCCAGACCTGTCTCTTATACACATCTCCGAGCCCACGAGACGATGCAGTATCTCGTATGCCGTCTTCGGCGTGAAAAGGGGGGG	FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFF,FFFFF:F,FF,F:,FF,F,F,,::FFF	NM:i:0	MD:Z:75	MC:Z:77S74M	AS:i:75	XS:i:69	XA:Z:chr5,-181470707,76S69M6S,0;chr2,+113615216,6S69M76S,0;chr1,+349618,6S69M76S,0;
A00928:207:HYLCHDSXY:2:2602:32705:31892	99	chr5	15259	60	151M	=	15692	584	GTTCACAGCCTTAAAGGTTAGTAGAGGTTTACCACACAGGAGGTAGAGGTTGCAGTGAGCTGAGATTGCACCTTTGCACTCCAGCCTGGGCGACAGAGCCAGACCCTGTCTCAAAAAAAAATTTTTTTAAAGGAAAACTATAGCTATTGTG	F,FFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFF,FFFFFFF,FFF,FFFF	NM:i:0	MD:Z:151	MC:Z:151M	AS:i:151	XS:i:102

比如这个是我们之前的比对的情况。从flag标记中我们可以看到，有的是第一条链上的比对到了正链，有的是第二条链上的比对到了正链。我们如何去解释这种情况呢？
而且这种情况还并不少见。我想和师姐交流一下，稍后。

6。对于出现在负链上的repeatfamily的家族的坐标，我们如何转换成正链的坐标进行处理？

如之前#4中提到的。对于比对到负链上的reads，在sam文件的结果中，我们将其转换为了在正链上的互补链的序列的信息。故而对于坐标而言，我想也是这条互补链与参考基因组比对的时候的最左边的坐标。也就是说在这里，所有的坐标都转换为了与参考基因组相对的坐标。
但在咱们repeatfamily的序列的文件中，显示的负链的坐标是什么呢？（是从5’-3’的负链的坐标？还是3’-5’的与正链相对的坐标？）

#bin	swScore	milliDiv	milliDel	milliIns	genoName	genoStart	genoEnd	genoLeft	strand	repName	repClass	repFamily	repStart	repEnd	repLeft	id
0	1892	83	59	14	chr1	67108753	67109046	-181847376	+	L1P5	LINE	L1	5301	5607	-544	1
1	2582	27	0	23	chr1	8388315	8388618	-240567804	-	AluY	SINE	Alu	-15	296	1	1
1	4085	171	77	36	chr1	25165803	25166380	-223790042	+	L1MB5	LINE	L1	5567	6174	0	4
1	2285	91	0	13	chr1	33554185	33554483	-215401939	-	AluSc	SINE	Alu	-6	303	10	6
1	2451	64	3	26	chr1	41942894	41943205	-207013217	-	AluY	SINE	Alu	-7	304	1	8
1	1587	272	100	49	chr1	50331336	50332274	-198624148	+	HAL1	LINE	L1	773	1763	-744	9
1	1393	280	82	51	chr1	58719764	58720546	-190235876	+	L2a	LINE	L2	2582	3418	-8	1
2	5372	165	14	27	chr1	75496057	75497775	-173458647	+	L1MA9	LINE	L1	5168	6868	-30	1
2	536	349	146	56	chr1	92274205	92275925	-156680497	+	L2	LINE	L2	406	2306	-1113	1

所以，我们一般情况下，这里的负链上的坐标指的是什么？

这里有一个细节要注意：
（1）这边的genoStart是0-based，也就是从0开始计数。0,1,2,3，……
（2）而sam文件中的坐标是1-based，也就是从1开始计数。1,2,3,4,5……
所以这也是为什么会有一位的差别。

从这张图上看，我好像明白了。这里的序列的genoStart和genoEnd好像是和参考基因组相比较的那个坐标。而且从文本中看，这些family好像是按照坐标的信息从小到大进行排列的。
如果果真是这样，那么这个事情就会变得简单很多。把负链的family按照正链同样处理即可。可能在最后的结果中，不同的family中还有一些重合的reads。就是这条reads既落在这个family的范围内，又落在那个family的范围内。因为我们本身是一个近似的结果。
这个时候还需要考虑比对的正负链了吗？
就是说是计算比对到负链上reads？
我还是觉得不用。因为我们这次计数的目的就是统计那块区域染色体是开放的。如果染色体的这块区域，可以检测到其正链的reads，那势必负链也是开放的（这个有依据吗？）。
所以，我们仅仅需要统计正链和负链比对到染色体的那块横跨的范围就可以了？

好像还是有点问题，需要再和师姐讨论明确一下：
（1）R1可能同时检测到既可以比对到正链，也可以比对到负链的reads嘛？
（2）ATAC-seq数据只关注染色体开放的区域，对于开放的正负链会有可能观测到嘛？因为转录的过程，的确是在一条链上进行的。但是咱们不是有两条同源染色体嘛？可能同时是不同的编码链转录嘛？

==>这个问题，事关我们如何设计我们的统计代码？

7。思考如何统计不同的细胞类型中的family的问题。

既然是不同的细胞类型，就势必会需要将这所有的细胞的reads先拆分成单个的细胞。但是，目前面临的问题可能是：
（1）首先我不知道如何拆分成单个细胞：无论是根据barcode还是barcode+spacer，得到的结果都和我们的预期10000左右细胞并不相同。另外，我也不太明白index这个序列的有什么用处（在scRNA中的index是用来标记样本的，比如来自不同的人）。但是这里的barcode和index都挺让我不解的。barcode相同，index不同。index相同，barcode不同。

===>关于不理解的细节可以再继续的整理一下。最重要的还是和师姐多多讨论这件事。

（2）另外可能需要问师兄要一下这方面的细胞类型的注释信息。
（3）我将来自于相同细胞类型的reads整合在一起，然后和参考基因组进行比对（或者是和subfamily家族的参考序列），通过类似的方式，看某一特定类型的细胞，在我们感兴趣的region，是否开放？另外，就是不同的family在不同的细胞中的表达的偏好性。我们或许可以找到一些特异调控的family，想想都觉得很有意思。
（4）计算量好大啊！到目前为止五号染色体的region还没有跑完。而我们就相当于不同细胞类型的result.sam都要进行这样的一个过程。让人忧伤。
（5）也担心会发现很少，因为从整体看就不多。

8。比较使用不同的方法得到的我们感兴趣的区域的reads的条数是否相同？如果不同，为什么？

写这个的过程，其实是对我们的思路进行整理的过程。
方案一： 将我们的ATAC-seq的数据和重复序列的数据库进行比较。将比对到的和我们感兴趣那一段的reads相同的家族的序列再次提取出来（这个时候可能就掩盖了双端序列的这一层信息），在与参考基因组进行匹配。然后再从中提取匹配到5号染色体上，我们感兴趣的那个region的reads的名字。
方案二： 将我们的ATAC-seq的数据与参考基因组进行比较。首先将双端最优匹配到五号染色体上的reads取出来。然后从这匹配到5号染色体的reads中，根据坐标提取，比对到我们感兴趣的区域的reads的名字。

我们想一下，这两种方法得到的结果之间是否有overlap的地方。目前方案二，我是得到了比较完整的结果。现在就是方案一的结果的提取，遇到了一些问题。我现在继续去处理这件事。

awk '$4>80029677 && $4<80029976' AluYe2_chr5.sam >AluYe2_chr5_interest.sam

A00928:207:HYLCHDSXY:2:1144:25898:22075	0	chr5	80029767	60	151M	*	0	0	GGCTAACATGGTGAAACCCCGTCTCTACTAAAAATACAAAAAAAAAAAAAAAAATTAGCCGGGCATGGTGGTGGGCACCTTTAGTCCCAGCTACTCGGGAGGCTGAGATAGGAGAATGGCGTGAACCCGGGAGGCGGAGGTTGCAGTAAGC	*	NM:i:0	MD:Z:151	AS:i:151	XS:i:97
A00928:207:HYLCHDSXY:2:2136:5195:35603	0	chr5	80029771	6	151M	*	0	0	AACATGGTGAAACCCCGTCTCTACTAAAAATACACACACACAAAAAAAAATTAGCCGGGCATGGTGGTGGGCACCTTTAGTCCCAGCTACTCGGGAGGCTGAGATAGGAGAATGGCGTGAACCCGGGAGGCGGAGGTTGCAGTGAGCCGAG	*	NM:i:5	MD:Z:34A1A1A1A102A7	AS:i:126	XS:i:121	XA:Z:chr8,-47957577,110M7I34M,9;
A00928:207:HYLCHDSXY:2:2136:5195:35603	0	chr5	80029810	25	151M	*	0	0	ACAAAAAAAAATTAGCCGGGCATGGTGGTGGGCACCTTTAGTCCCAGCTACTCGGGAGGCTGAGATAGGAGAATGGCGTGAACCCGGGAGGCGGAGGTTGCAGTGAGCCGAGATCGCGCTGCTGCACTCCAGCCTGGGCGACAGCGAGACT	*	NM:i:2	MD:Z:1A102A46	AS:i:144	XS:i:132	XA:Z:chr8,-47957536,7M2D142M2S,4;
A00928:207:HYLCHDSXY:2:1144:25898:22075	0	chr5	80029821	50	151M	*	0	0	TTAGCCGGGCATGGTGGTGGGCACCTTTAGTCCCAGCTACTCGGGAGGCTGAGATAGGAGAATGGCGTGAACCCGGGAGGCGGAGGTTGCAGTAAGCCGAGATCGCGCTGCTGCACTCCAGCCTGGGCGACAGCGAGACTCCATCTCAAAA	*	NM:i:0	MD:Z:151	AS:i:151	XS:i:128	XA:Z:chr8,-47957525,18M2D133M,5;

忽略掉双端的情况，也就是说在我们重复序列的比对的结果中，有两条reads是比对到了我们感兴趣的区域中。

A00928:207:HYLCHDSXY:2:1144:25898:22075
A00928:207:HYLCHDSXY:2:2136:5195:35603

然后我们接着检验这两条reads在不在我们之前从染色体层面的比较的结果中。其次这两条reads在序列上合不合我们感兴趣的区域有overlap的情况（也就是说，是否真的比对上了）。

能，感动哭了。


也就是说我们的这两种方法之间是可以相互印证的。
那我想知道，我们现在通过坐标找到的序列是否在原先的重复序列的list中呢？也想做一个交叉验证。

(base) [xxzhang@mu02 trash]$ grep "A00928:207:HYLCHDSXY:2:2216:1931:8437" AluYe2_chr.sam                                                                                                                                                           
(base) [xxzhang@mu02 trash]$ grep "A00928:207:HYLCHDSXY:2:2259:28085:19319" AluYe2_chr.sam                                                                                                                                                            
(base) [xxzhang@mu02 trash]$ grep "A00928:207:HYLCHDSXY:2:2419:4399:22200" AluYe2_chr.sam

那么，如果这部分的序列没有比对到AluYe2这个家族上，又会比对到哪里呢？

(base) [xxzhang@mu02 family]$ grep "A00928:207:HYLCHDSXY:2:2216:1931:8437" redup.F2.sam
A00928:207:HYLCHDSXY:2:2216:1931:8437   99      AluYh9  1       0       70S81M =1       96      TTACTAATCGTCAGTGCCTAATTGCTAAAGCAAAATATCTTCTAAATTCTGGGGAGGAAATTCTACTTTGGGCCAGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGGGAGGCCGAGGCGGGCTGATCGTGAGGTCAGGAGATCGAG FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF::FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFAS:i:134 XS:i:134        XN:i:0  XM:i:4  XO:i:0  XG:i:0  NM:i:4  MD:Z:4G42A9G4A18YS:i:164        YT:Z:CP MQ:i:0  MC:Z:55S96M     ms:i:5416
A00928:207:HYLCHDSXY:2:2216:1931:8437   147     AluYh9  1       0       55S96M =1       -96     GCCTAATTGCTAAAGCAAAATATCTTCTAAATTCTGGGGAGGAAATTCTACTTTGGGCCAGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGGGAGGCCGAGGCGGGCTGATCGTGAGGTCAGGAGATCGAGACCATCCTGGCTAAC FFFF:FFFFFFF,FFFFFFFFF:FFFF,FFFFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF::FFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFAS:i:164 XS:i:164        XN:i:0  XM:i:4  XO:i:0  XG:i:0  NM:i:4  MD:Z:4G42A9G4A33YS:i:134        YT:Z:CP MQ:i:0  MC:Z:70S81M     ms:i:5551
(base) [xxzhang@mu02 family]$ grep "A00928:207:HYLCHDSXY:2:2259:28085:19319" redup.F2.sam
(base) [xxzhang@mu02 family]$ grep "A00928:207:HYLCHDSXY:2:2419:4399:22200" redup.F2.sam
#有些甚至都没有比对上

(base) [xxzhang@mu02 family]$ grep "A00928:207:HYLCHDSXY:2:2419:4399:22200" result_HumSub.sam
A00928:207:HYLCHDSXY:2:2419:4399:22200  81      AluYe2  209     2       87S64M  =       210     0       CCCCCCCCCCCCCCCCCCTTTTCAAGCAGAAGACGGCATACGAGATGAAAGTCTGTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGAGGCGGAGGTTGCAGTAAGCCGAGATCGCGCTGCTGCACTCCAGCCTGGGCGACAGCGAGACT      FFFFFFFFFFFFFF:F,,:,:FFF:FF,::FFFFFFFFFFFF::FFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF AS:i:100        XN:i:0  XM:i:4  XO:i:0  XG:i:0  NM:i:4  MD:Z:9C7G14C0A30     YT:Z:UP
A00928:207:HYLCHDSXY:2:2419:4399:22200  161     AluYe2  210     2       70M81S  =       209     0       AGGCGGAGGTTGCAGTAAGCCGAGATCGCGCTGCTGCACTCCAGCCTGGGCGACAGCGAGACTCTGTCTCTTATACACATCTGACGCTGCCGACGACAGACGCGTAAGGTTACCTAGCGAGTGTAGATCTCGGTGGTCGCCGTATCATTAA      FFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFF,,FFF:,FF AS:i:105        XN:i:0  XM:i:5  XO:i:0  XG:i:0  NM:i:5  MD:Z:8C7G14C0A31C5   YT:Z:UP

(base) [xxzhang@mu02 family]$ grep "A00928:207:HYLCHDSXY:2:2259:28085:19319" result_HumSub.sam
A00928:207:HYLCHDSXY:2:2259:28085:19319 99      AluYb8a1        187     0       4S53M94S        =       187     -240    AGATAGGAGAATGGCGTGAACCCGGGAGGCGGAGGTTGCAGTAAGCCGAGATCGCGCTGCTGCTGTCTCTTATACACATCTCCGAGCCCACGAGACAGACTTTCATCTCGTATGCCGTCTGCTGCTTGAAAGGGGGGGGGGGGGGGGGGGG      :FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFF:,FFF,FFFF:FF::FFFF,:F,,F,,F,,FFFF:FF:FFFFFFF AS:i:92 XS:i:92 XN:i:0  XM:i:2  XO:i:0  XG:i:0       NM:i:2  MD:Z:30C7G14    YS:i:92 YT:Z:CP
A00928:207:HYLCHDSXY:2:2259:28085:19319 147     AluYb8a1        187     0       93S53M5S        =       187     240     TTTAATGATACGGCGACCACCGAGATCTACACACAACGGTCATAGCTTCGCGTCTGTCGTCGGCAGCGTCAGATGTGTATAAGAGACAGAGATAGGAGAATGGCGTGAACCCGGGAGGCGGAGGTTGCAGTAAGCCGAGATCGCGCTGCTG      FF:FFFFF,,FFFFFFFFFFFFFFF:FFFF:F:FFFFFFF:F:FFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF AS:i:92 XS:i:92 XN:i:0  XM:i:2  XO:i:0  XG:i:0       NM:i:2  MD:Z:30C7G14    YS:i:92 YT:Z:CP

总结一下，没有对上的情况是这样的。要么是这些reads在序列上匹配到了其他的家族（因为毕竟都属于一个subfamily，在序列上具有很大的相似性）。要么是比对到了这个家族，但是在我一系列的筛选的过程中，被筛去了。故而，造成了reads数量的减少。

===>因此现在比较起来看，老师的想法还是很有道理的，结果也的确证实为这样的。
总结一下：subfamily方法得到的reads一定在我们用坐标得到的结果中。而用坐标方法得到的reads不一定在这个subfamily的结果中，因为可能比对的过程中被筛去，或者因为序列的极其相似，所以比对到了其他的家族中。

这个问题也论证结束了。