实验记录 | scATAC-seq数据的比对（三）

1。subfamily计数

接着咱们之前的一些处理的步骤。就是说写一个指令，遍历执行语句：

samtools view redup.F2.bam | grep "AluYe2" >AluYe2.sam
cat AluYe2.sam | wc -l
#获取屏幕的输出的结果，格式为：AluYe2\tnumber(到时候直接画成histgram)
#甚至觉得有了那个list之后，可以直接用代码搞定。
#这里要替换的值，可能只有subfamily的名字（如果让我一个个手动去看的话，鬼知道要猴年马月？）

先使用R提取了那个subfamily_name的list。然后使用perl批量的执行。
后来就写了个代码，批量的循环解决了（就是懒啊）。我爱写代码的一个重要的原因，就是它可以让我远离重复性的劳动。而且即使再细致的人工，始终都没有计算机精确。这也是我迷信它的原因。只要方法正确，就一定可以达到我们想要的结果。

#!perl
open (MARK, "< subfamily_name.txt") or die "can not open it!";
while ($line = <MARK>){ 
        $line=~s/\r\n//g;
        system_call("samtools view redup.F2.bam | grep "."\"".$line."\""." >".$line.".sam");
	print $line,":";
	system_call("cat ".$line.".sam  | wc -l");
		
#获取屏幕的输出的结果，格式为：AluYe2\tnumber(到时候直接画成histgram)
#甚至觉得有了那个list之后，可以直接用代码搞定。
} 
#print @ID;
close(MARK);

sub system_call
{
  my $command=$_[0];
  print "\n\n".$command."\n\n";
  system($command);
}

最后将屏显的结果，放到notepad上进行初步的筛选。
（1）首先使用正则表达式^.*cat.*\r?\n和^.*samtools.*\r?\n删除有这个的行。
（2）然后将换行符\r\n替换为空。
经过上述步骤的筛选，得到了我们想要的结果。

最后就是将这个list提交到R中，绘制barplot。
（反正就是，一开始想哪一个语言能够解决这个事儿，不影响彼此之间交互替换着使用）

2。与参考基因组进行比对

 bwa mem ./hg38.fasta aluY.fasta >aluY_chr.fasta
 mv aluY_chr.fasta aluY_chr.sam
 grep "hg38" aluY_chr.sam

@PG ID:bwa PN:bwa VN:0.7.17-r1188 CL:bwa mem ./hg38.fasta aluY.fasta
hg38_rmsk_AluY 0 chr5 80029678 60 299M * 0 0 GGCCAGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGGGAGGCCGAGGCGGGCTGATCGTGAGGTCAGGAGATTGAGACCATCCTGGCTAACATGGTGAAACCCCGTCTCTACTAAAAATACAAAAAAAAAAAAAAAAATTAGCCGGGCATGGTGGTGGGCACCTTTAGTCCCAGCTACTCGGGAGGCTGAGATAGGAGAATGGCGTGAACCCGGGAGGCGGAGGTTGCAGTAAGCCGAGATCGCGCTGCTGCACTCCAGCCTGGGCGACAGCGAGACTCCATCTCAAAAAAAAA * NM:i:0 MD:Z:299 AS:i:299 XS:i:245 XA:Z:chr8,-47957520,23M2D133M7I136M,15;

最终是完全匹配到了第5号染色体上。
我现在想做一件有意思的事情。
在bowite比对的时候，我们的reads不是比对到了aluYe2这个subfamily上嘛。我现在就想看看，把这部分的subfamily提取出来，确定他们染色体的位置。
先提取sam文件（之前计数的时候已经提取过了），然后可能需要一个sam转成fastq的小工具。

遇到个问题，原来我们直接检索的那个sam文件没有SQ的头。
解决思路来源链接：https://blog.csdn.net/g_r_c/article/details/8252256

 samtools faidx humsub.fasta
 samtools view -bt humsub.fasta.fai AluYe2.sam AluYe2.bam

即完成转换。
然后将bam转换为fastq文件。
解决思路来源链接：https://www.jianshu.com/p/1a3e090805cd
（由此可见awk真的好有用啊，我要学会它！）

samtools view AluYe2.bam -O SAM | awk -F"\t" '{print ">"$1"\n"$10}' > AluYe2.fa
head AluYe2.fa

屏幕输出：

A00928:207:HYLCHDSXY:2:1101:1217:17002
ACATAGTTGGAAGTAAAGCACTCCTCAGCAAATGTAAAAGAACAGAAATCACAACAAACTGTCTCTCAGACCACAGTGCAATCAAATTAGAACTCAGGATTAAGAAACTCACTCAGGGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGC
A00928:207:HYLCHDSXY:2:1101:12183:3709
TGCCCACAGGACTGAATCCATTACAGGTATAAAATATAAAAAGTGCTTTAAATTTTTCGAAACCTGGCTGGGCACGGTGGCTCACGCCTGTAATCCCAGCACTTTGGGAAGCTGAGGCACTGTCTCTTATACACATCTCCGAGCCCACGAG
A00928:207:HYLCHDSXY:2:1101:15682:29058
GACCAATGACTTAACGACTCAGATTAAAGGGCTGGAATATCTCATTCCTATCCCTTGCAAGATCCCTGAGGAGATCAGCAAACTATTAAAGAAATTTAAGGCTGGGCGTGGTGGCTCATGCCTGTAATCCCAGCACTTTGGGAGGCCGAGG
A00928:207:HYLCHDSXY:2:1101:18810:35196
GCCTCAAGTAACATAATACAAAAATAGTCAAAATTTTCCTGCTTTGAAACAATCTACTTAATTGCATGCACAGTTTTCTGAGTCTGCTTTGAAATTAAAAATAAAGTGTGGGCCGGGCACGGTGGCTCACGCCTGTAATCCCAGCACTTTG
A00928:207:HYLCHDSXY:2:1101:23104:15702

然后我们拿这部分数据和我们的参考基因组进行比对。

 bwa mem ./hg38.fasta AluYe2.fa >AluYe2_chr.sam
 grep "chr5" AluYe2_chr.sam |wc -l
3142
cat AluYe2_chr.sam | wc -l
54452 #把运动锻炼当做是一种生活。